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A 


ABSTRACT 


A  mathematical  model  of  a  multicomponent  reacting  nonsimilar  laminar  or  turbulent 
boundary  layer  flow  including  transverse  curvature  effects  is  presented,  and  a 
method  of  solution  is  given.  The  formulation  is  an  extension  of  an  earlier  work 
in  which  thermal  and  molecular  diffusion  were  treated  in  terms  of  a  bifurcation 
approximation  for  binary  diffusion  coefficients.  In  the  present  analysis,  a  tur¬ 
bulent  model  is  added  which  employs  a  mixing  length  model  for  eddy  viscosity  in 
the  wall  region  with  consideration  of  injection  or  suction  effects.  The  wake  re¬ 
gion  eddy  viscosity  is  taken  to  be  proportional  to  the  free  stream  velocity  and 
local  velocity  defect  thickness.  Transverse  curvature  effects  are  also  incorpo¬ 
rated  into  the  present  analysis.  A  modification  of  the  Levy-Lees  transformation 
is  used  to  transform  the  equations  of  motion  to  the  (£,n)  coordinate  plane,  where 
the  conservation  equations  are  integrated  across  boundary  layer  strips.  Deriva¬ 
tives  in  the  normal  direction  are  related  to  one  another  by  Taylor  series  trun¬ 
cated  to  reflect  a  quadratic  or  cubic  approximation,  and  streamwise  derivatives 
are  expressed  in  finite  difference  form.  The  resultant  set  of  equations  is  solved 
by  general  Newton-Raphson  iteration. 

(Distribution  Limitation  Statement  No.  2) 
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SYMBOLS 


parameter  used  in  the  solution  of  the  mixing  length  equation  (de¬ 
fined  by  Equation  (131)) 


parameter  used  in  the  solution  of  the  mixing  length  equation  (de¬ 
fined  by  Equation  (132)) 


constant  introduced  in  the  a„  constraint  (Equation  (58)) 

n 

constant  introduced  in  the  approximation  for  multicomponent  ther¬ 
mal  diffusion  coefficients  embodied  in  Equation  (25) .  Tentatively 
established  by  correlation  of  data  to  be  -0.5 


product  of  density  and  viscosity  normalized  by  their  reference 
values  (defined  by  Equation  (66)) 


frozen  specific  heat  of  the  gas  mixture  (defined  by  Equation  (17)) 


property  of  the  gas  mixture  which  reduces  to  Cp  when  diffusion  co¬ 
efficients  are  assumed  equal  for  all  species  L  (defined  by  Equa¬ 
tion  (26) ) 


specific  heat  of  species  i 


coefficients  defined  in  finite-difference  representation  of  stream- 
wise  derivatives  (defined  in  Equations  (112)  and  (113)  for  two- 
and  three-point  difference  relations,  respectively) 


a  reference  binary  diffusion  coefficient  introduced  by  the  approx¬ 
imation  for  binary  diffusion  coefficients  embodied  in  Equation  (19) 


multicomponent  thermal  diffusion  coefficient  for  species  i 


multicomponent  diffusion  coefficient  for  species  i  and  j 


diffusion  coefficient  for  all  species  when  all  A  ^  are  equal 
binary  diffusion  coefficient  for  species  i  and  j 


errors  for  the  various  equations  during  Newton-Raphson  iteration 
(driven  toward  zero  in  the  iteration) 


stream  function  (defined  by  Equation  (59)) 


diffusion  factor  for  species  i  introduced  by  the  approximation  for 
binary  diffusion  coefficients  embodied  in  Equation  (19) 
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SYMBOLS 

(continued) 


h  static  enthalpy  of  the  gas  (defined  by  Equation  (15)) 

hw  static  enthalpy  of  the  gas  at  the  wall 


property  of  the  gas  mixture  which  reduces  to  the  static  enthalpy 
h  when  diffusion  coefficients  are  assumed  equal  for  all  species 
(defined  by  Equation  (27)) 


enthalpy  of  surface  material  (a.g.,  char)  removed  by  combustion, 
sublimation,  or  vaporization 


enthalpy  of  gas  which  enters  boundary  layer  without  phase  change 
at  the  surface,  (e.g.,  pyrolysis  gases) 


enthalpy  of  species  i  (defined  by  Equation  (16)) 


heat  of  formation 


enthalpy  of  component  surface  material  (e.g.,  silica)  removed 
in  the  condensed  phase  (e.g.,  by  melting  with  subsequent  liquid 
runoff  or  by  spallation) 


total  enthalpy  (defined  by  Equation  (14)) 


diffusional  mass  flux  of  species  i  per  unit  area  away  from  the 
surface 


diffusional  mass  flux  of  element  k  per  unit  area  away  from  the 
surface 


total  number  of  elements?  also  mixing  length  constant 


mass  fraction  of  molecular  species  i 


k 


total  mass  fraction  of  element  (or  base  gas)  k  contained  in  sur¬ 
face  material  (e.g.,  char)  removed  by  combustion,  sublimation,  or 
vaporization 


k 


total  mass  fraction  of  element  (or  base  gas)  k  contained  in  gas 
which  enters  boundary  layer  without  phase  change  at  the  surface 
(e.g.,  pyrolysis  gases) 


total  mass  fraction  of  element  (or  base  gas)  k  irrespective  of  mo¬ 
lecular  configuration  (defined  by  Equation  (lr‘)) 
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SYMBOLS 

(continued) 


£  mixing  length  (defined  by  Equation  (55)) 


£  dimensionless  mixing  length  (defined  by  Equation  (71)) 


L  parameter  used  in  mixing  length  formulation  (defined  by  Equation 

(128)) 


m  mass  flow  rate  per  unit  area 

mc  mass  removal  rate  per  unit  area  of  surface  material  (e.g.,  char) 

by  combustion,  sublimation,  or  vaporization 

m  mass  flow  rate  per  unit  area  of  gas  which  enters  boundary  layer 

g  without  phase  change  at  the  surface  (e.g.,  pyrolysis  gases) 

*  til 

mr  mass  removal  rate  per  unit  area  of  £  component  surface  material 

r£  (e.g.,  silica)  in  the  condensed  phase  (e.g.,  by  melting  with  sub¬ 

sequent  liquid  runoff  or  by  spallation) 


%  molecular  weight  of  the  gas  mixture 

U\^  molecular  weight  of  species  i 

N  number  of  nodal  points  across  the  boundary  layer  selected  for  the 

purpose  of  the  numerical  solution  procedure 

p  dummy  variable  representing  f',  H^,  or  Kk 

P  pressure;  also  a  parameter  used  in  the  mixing  length  formulation 

(defined  by  Equation  (125)) 


partial  pressure  of  species  i 

Pr  frozen  Frandtl  number  of  the  gas  mixture  (defined  by  Equation  (86)) 


Prt  turbulent  Prandtl  number  (defined  by  Equation  (69)) 

q  diffusional  heat  flux  per  unit  area  away  from  the  surface 

£1 


gcond 


heat  conduction  per  unit  area  into  the  surface  material 


one-dimensional  radiant  heat  flux  (toward  the  surface),  than  is, 
the  net  rate  per  unit  area  at  which  radiant  energy  is  transferred 
across  a  plane  in  the  boundary  layer  parallel  to  the  surface 
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SYMBOLS 

(continued) 


local  radius  in  the  boundary  layer  in  a  meridan  pl'ne  for  an  axi- 
symmetric  shape 

local  radius  of  body  in  a  meridian  plane  for  an  axisymmetrie  shape 
universal  gas  constant 

Reynolds  number;  subscripted  with  the  length  scale  if  other  than  s 

effective  nose  radius  for  Newtonian  flow 

distance  along  body  from  stagnation  point  or  leading  edge 

reference  system  Schmidt  number  (defined  by  Equation  (88)) 

turbulent  Schmidt  number  (defined  by  Equation  (68)) 

parameter  defined  to  simplify  problems  with  transverse  curvature; 
see  Equation  (55) 

static  temperature 

velocity  component  parallel  to  body  surface 
shear  velocity,  defined  in  Equation  (49) 
velocity  component  normal  to  body  surface 
mole  fraction  of  species  i 

truncated  series  obtained  in  Taylor  series  expansion  of 


J1  f'p  dn  (defined  by  Equation  (117)) 
i-1 


distance  from  surface  into  the  boundary  layer,  measured  normal  to 
the  surface 


dimensionless  y-coordinate  defined  by  Equation  (49)) 


constant  in  the  mixing  length  differential  equation  (see  Equation 
(45)) 
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SYMBOLS 

(continued) 


Z.  a  quantity  for  species  i  which  is  introduced  as  a  result  of  the 

~  approximation  for  binary  diffusion  coefficients  and  reduces  to 

when  all  diffusion  coefficients  are  assumed  equal  (defined  by 
Equation  (20)) 


Z.  a  quantity  for  element  (or  base  species)  k  v/hich  is  introduced  as 

*  a  result  of  the  approximation  for  binary  diffusion  coefficients 

and  reduces  to  when  all  diffusion  coefficients  are  assumed 
equal  (defined  by  Equation  (28)) 


ZP^,ZP2/ • • • 


a* 


aki 


«,A&-1 


Af^Aq,... 


truncated  series  obtained  in  Taylor  series  expansion  of  integrals 
involving  nonsimilar  terms  (defined  by  Equation  (124)) 

flux  normalizing  parameter  (defined  by  Equation  (82)) 

normalizing  parameter  used  in  definition  of  n  (see  Equation  (59)) 
defined  implicitly  by  use  of  a  constraint  such  as  Equation  (60) 

mass  fraction  of  element  (or  base  species)  k  in  species  i 

streamwise  pressure-gradient  parameter  (defined  by  Equation  (67)) 

y-dimension  normalizing  parameter  (defined  by  Equation  (70)) 

logarithmic  distance  between  two  streamwise  positions  denoted  by 
the  subscripts  Z  and  2,-1  (defined  by  Equation  (114)) 

corrections  for  q,  q,...,  during  Newton-Raphson  iteration 


6*  displacement  thickness  (defined  by  Equation  (52)) 

<5*  incompressible  or  velocity  displacement  thickness  (defined  by 

1  Equation  (53)) 

6ri  distance  between  two  boundary  layer  nodal  points  (defined  by 

Equation  (105) ) 

n,g  transformed  coordinate  in  a  direction  normal  to  the  surface  (de¬ 

fined  by  Equation  (61)).  Note:  the  hat  is  dropped  from  n  through¬ 
out  most  of  the  report 


6  angle  between  a  surface  normal  and  a  normal  to  the  body  center- 

line;  also  time  in  discussions  of  a  charring  ablation  program 


\ 


thermal  conductivity 
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SYMBOLS 

(continued) 


shear  viscosity 


properties  of  the  gas  mixture  (defined  by  Equations  (21)  and  (24)) 
which  reduce  to  unity,  to  %  to  1/2?;,  and  to  In  %  respectively, 
for  assumed  equal  diffusion  coefficients 


kinematic  viscosity 


transformed  streamwise  coordinate  (defined  by  Equation  (61)).  Note 
the  hat  is  dropped  from  t,  throughout  most  of  the  report 


density 


total  mass  flux  per  unit  aiaa  into  the  boundary  layer 


individual  species  turbulent  eddy  diffcsivity  (defined  by  Equation 
(2)) 


average  turbulent  eddy  diffusivity,  where  it  is  assumed  that  all 
peDi  =  P£D 

turbulent  eddy  conductivity  (defined  by  Equation  (18)) 


turbulent  eddy  viscosity  (defined  by  Equation  (12) 


dimensionless  eddy  viscosity  (defined  by  Equation  (73)) 


Stefan-Boltzmann  constant 


local  shear  stress 


"elemental" source  term  (defined  by  Equation  (29)) 


rate  of  mass  generation  of  species  i  per  unit  volume  due  to  chemi¬ 
cal  reaction 


Subscripts 


pertains  to  boundary-layer  edge 


equil 


pertains  to  surface  equilibrium  requirement 
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SYMBOLS 

(concluded) 


i  pertains  to  the  species  or  to  the  ifc^  nodal  point  in  the  bound¬ 

ary  layer,  starting  with  i  =  1  at  the  surface 

j  pertains  to  jth  species 

til 

k  pertains  to  k1"  element  (or  base  species) 

i  j. 

l  pertains  to  *  streamwise  position 

th 

m  pertains  to  m  iteration  during  the  Newton-Raphson  iteration 

process 

n  pertains  to  the  n*"*1  nodal  points ,  corresponding  to  the  outer  edge 

of  the  boundary  layer  solution 

sp  pertains  to  the  stagnation  point 

s.s.  pertains  to  the  steady  state  energy  balance  requirement 

w  pertains  to  wall 

1  reference  condition,  usually  taken  as  zero  streamline  from  invis- 

cid  solution  (synonymous  with  boundary- layer  edge  in  the  absence 
of  an  entropy  layer) 

Superscripts 

k  equal  to  unity  for  axisymmetric  bodies  and  zero  for  two-dimensional 

bodies 


signifies  that  quantity  is  normalized  by  a*  (e.g.,  j£  =  j^/ct*) 

Represents  partial  differentiation  with  respect  to  n  or  n  (usually 
ri  unless  otherwise  noted) .  Represents  turbulent  fluctuation  in 
Section  II. 
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SECTION  I 
INTRODUCTION 

A  computational  procedure  is  described  which  is  suitable  for  obtaining  accu¬ 
rate  numerical  solutions  of  the  nonsimilar  multicomponent  laminar  and/or  turbulent 
boundary  layer  with  arbitrary  equilibrium  or  nonequilibrium  chemical  systems,  un¬ 
equal  diffusion  and  thermal  diffusion  coefficients  for  all  species,  radiation 
absorption  and  emission,  second  order  transverse  curvature  effects,  and  a  variety 
of  surface  boundary  conditions  including  intimate  coupling  with  transient  char¬ 
ring-ablation  energy  and  mass  balances.  A  Fortran  IV  computer  program  has  been 
developed  in  accordance  with  this  analysis  with  the  exceptions  that  1)  the  chemi¬ 
cal  system  is  presently  limited  to  equilibrium  in  the  boundary  layer,  with  or 
without  selected  rate-controlled  surface  reactions  or  surface  catalyzed  reactions, 
and  2)  radiation  absorption  and  emission  within  the  boundary  layer  is  not  permit¬ 
ted  in  the  version  of  the  program  reported  here.  This  computer  program,  desig¬ 
nated  BLIMP,  for  Boundary  Layer  Integral  Matrix  Procedure,  is  described  in  ref¬ 
erence  1.  The  analysis  and  computer  program  described  herein  are  extensions  of 
the  previously  developed  BLIMP  program  for  laminar  boundary  layer  flows  described 
in  reference  2.  The  turbulent  model  which  has  been  incorporated  was  reported 
earlier  with  a  restriction  to  incompressible  flows  in  reference  3. 

The  computational  procedure  was  developed  while  attempting  to  take  advantage 
of  the  most  attractive  features  of  other  boundary- layer  procedures.  In  light  of 
the  application  of  the  procedure  to  be  adopted,  certain  specific  requirements 
were  appropriate.  In  particular,  minimization  of  the  number  of  "nodal  points" 
required  to  obtain  a  solution  was  judged  to  be  of  prime  importance  as  a  conse¬ 
quence  of  the  relatively  large  times  associated  with  state  calculations  for  a  gen¬ 
eral  chemical  environment  and,  in  the  streamwise  direction,  because  of  the  desire 
to  couple  the  boundary  layer  procedure  to  a  transient  internal  conduction  or  abla¬ 
tion  solution. 

For  a  given  accuracy,  the  number  of  necessary  "nodal  points"  in  the  surface 
normal  direction  is  controlled  primarily  by  the  nature  of  the  functions  which  re¬ 
late  the  dependent  variables  (and  their  derivatives)  to  the  independent  variable. 
Thus  the  continuous  functions  typically  used  in  integral  relations  approaches  re¬ 
quire  fewer  "nodal  points"  than  the  functions  with  discontinuous  first  derivatives 
implied  by  most  finite  difference  approximations.  In  order  to  permit  relatively 
flexible  profiles,  sets  of  connected  quadratics  and  cubics  were  selected  to  rep¬ 
resent  enthalpy,  velocity,  and  elemental  concentrations.  The  first  derivatives 
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(and  second  derivatives  with  cubics)  of  these  functions  were  made  continuous  at 
the  connecting  points.  The  advantages  of  such  a  "spline  fit"  are  considered, 
for  example,  in  reference  4. 

If  the  general  integral  relations  approach  is  followed,  weighting  functions 
must  be  selected.  In  the  present  study,  as  in  reference  2,  step  weighting  func¬ 
tions  similar  to  those  used  by  Pallone  (reference  5)  were  used.  That  is,  the 
conservation  equations  are  integrated  between  nodal  points  (over  strips)  with  a 
unity  weighting  function. 

In  the  past  when  relatively  large  spacing  in  the  streamwise  direction  has 
been  desired,  iterative  procedures  have  generally  been  used  to  assure  accuracy 
and  stability.  Some  of  these  procedures  have  treated  the  solution  in  a  manner 
resembling  that  used  for  a  similar  solution  but  with  the  addition  of  finite  dif¬ 
ference  representations  for  the  nonsimilar  terms,  a  procedure  which  eliminates 
the  necessity  of  special  starting  techniques.  Using  this  basic  approach,  the 
specific  treatment  adopted  in  the  current  study  follows  most  closely  the  matrix 
procedure  used  by  Leigh  (reference  6)  wherein  the  iteration  is  a  consequence  of 
the  solution  of  a  set  of  linear  and  nonlinear  algebraic  relations .  The  general 
Newton-Raphson  technique  was  used  in  the  present  procedure  to  solve  these  simul¬ 
taneous  equations.  This  technique  results  in  linearized  coupling  between  all 
relations  required  to  characterize  the  boundary  layer,  and  thus  assures  a  more 
general,  rapid  and  stable  iterative  convergence. 

This  document  is  the  second  report  to  describe  completely  the  analysis  and 
solution  procedure  associated  with  the  BLIMP  program,  reference  2  being  the  first. 
The  addition  of  a  turbulent  boundary  layer  capability  and  transverse  curvature 
effects  to  the  program  provided  the  impetus  for  this  report,  therefore  these  cwo 
topics  will  receive  perhaps  a  disproportionate  share  of  discussion.  Much  of  the 
rest  of  the  analysis  and  solution  procedure  is  the  same  as  was  reported  in  ref¬ 
erence  2;  the  reader  is  referred  to  that  document  for  more  complete  discussions. 

This  report  concentrates  on  the  fluid  mechanical  aspects  of  the  problem 
and  describes  the  basic  numerical  solution  procedure.  The  procedures  employed 
for  calculating  the  equilibrium  state  of  the  gas  and  suggested  for  including 
rate-controlled  reactions  are  described  elsewhere  (reference  7)  since  they  are 
conveniently  treated  as  subroutines  to  the  basic  boundary  layer  computational 
procedure.  However,  the  terms  which  are  directly  involved  in  the  boundary  layer 
equations  such  as  the  "elemental  source  term"  which  arises  from  kinetic  consid¬ 
erations  are  included  in  the  present  development.  Similarly,  radiation  absorp¬ 
tion  and  emission  enters  directly  into  the  conservation  equations  only  as  a  net 
radiation  flux  term  in  the  energy  equation.  The  calculation  of  this  term  can  be 
conveniently  accomplished  by  a  subroutine.  Multicomponent  transport  properties 
are  based  on  the  approximation  reported  in  reference  8.  Modification  of  the 
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conservation  equations  as  a  consequence  of  this  approximation  is  described 
briefly  here,  and  more  completely  in  reference  2.  The  procedures  employed  for 
coupling  to  a  transient  charring  ablation  program  are  also  described  here. 

Section  II  includes  the  entire  mathematical  modeling  of  the  boundary  layer 
flow  including  discussions  of  the  general  conservation  equations,  turbulent  flow 
considerations,  transverse  curvature  effects,  coordinate  transformations,  and 
boundary  conditions.  Section  III  outlines  the  integral  matrix  method  for  solving 
the  simultaneous  differential  equations  including  the  integral  strip  relations, 
the  mixing  length  solution  procedure,  and  the  Newton-Raphson  iteration  technique. 
Section  IV  presents  some  of  the  results  obtained  with  the  program.  Section  V  de¬ 
scribes  the  CABLE  (Charring  Ablation  and  Boundary  Layer  Environment)  program 
wherein  the  boundary  layer  analysis  is  coupled  to  a  one-dimensional  charring 
ablation  analysis  at  each  body  station  and  presents  the  results  of  a  sample  run. 
Section  VI  contains  the  summary,  conclusions,  and  recommendations  for  further 
study. 
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SECTION  II 

MATHEMATICAL  MODEL  OF  THE  BOUNDARY  LAYER 


1.  GENERAL  CONSERVATION  EQUATIONS 

In  the  present  analysis,  the  usual  turbulent  flow  technique  of  breaking  the 
species,  velocity,  and  enthalpy  fields  into  mean  and  fluctuating  components,  time 
averaging,  and  making  appropriate  order  of  magnitude  approximations  is  used.  Tak¬ 
ing  the  results  of  these  manipulations  as  a  point  of  departure,  the  species  mass 
balance  equation  becomes 
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where  s  and  y  are  the  streamwise  and  normal  coordinates,  respectively,  u  and  v 
are  the  velocity  components  in  the  s  and  y  directions,  respectively,  is  the 
mass  fraction  of  species  i,  r  is  the  radius  from  the  body  centerline  to  the 
point  of  interest  in  a  meridian  plane  for  an  axisymmetric  shape,  <  is  zero  for  a 
flat  plate  and  unity  for  a  body  of  revolution,  p  is  the  density,  and  iJk  repre¬ 
sents  the  rate  of  mass  generation  of  species  i  per  unit  volume  due  to  chemical 
reaction.  The  individual  species  turbulent  eddy  diffusivity  peDi  is  defined  in 
terms  of  the  correlation  of  the  fluctuating  components  of  concentration  and  nor¬ 
mal  velocity,  that  is, 
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and  is  the  mass-diffusion  rate  of  species  i  due  to  molecular  processes.  Since 
transverse  curvature  is  to  be  included  in  the  present  analysis,  r  must  be  treated 
as  a  function  of  y  whereas  in  the  typical  boundary  layer  analysis,  r  is  set  equal 
to  rQ,  the  surface  radius.  The  relationship  between  r,  rQ,  and  y  is 

r(s,y)  =  r'Q(s)  +  y  cos  0  (3) 

Figure  1  helps  orient  the  reader  to  the  nomenclature  being  used. 

In  equation  (1)  and  in  other  conservation  equations  to  follow,  turbulent 
transport  terms  are  expressed  in  Boussinesq  form,  that  is,  eddy  viscosity,  eddy 
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diffusivity,  and  eddy  conductivity.  Hence  all  terms  are  time-averaged  quanti¬ 
ties  and  no  need  exists  for  using  a  superscript  bar.  In  the  order-of-magnitude 
arguments,  terms  of  the  following  types  have  been  eliminated:  1)  triple  corre¬ 
lations,  2)  derivatives  of  turbulent  correlations  parallel  to  the  wall,  and  3) 
correlations  involving  turbulent  components  of  molecular  transport  mechanisms. 

When  equation  (1)  is  summed  over  all  species,  the  global  continuity  equa¬ 
tion  results: 

3purK  ,  3pvrK 
3s  +  3y  0 

(4) 

Combining  equations  (1)  and  (4) ,  one  obtains  the  species  conservation  equation 


which  can  be  written  for  each  species  i  under  consideration.  The  molecular  dif¬ 
fusion  rate  is  expressed  in  general  as 


where  D^..  is  the  multicomponent  diffusion  coefficient  of  species  i  into  j,  dT 
is  the  multicomponent  thermal  diffusion  coefficient  of  species  i,  is  the  local 
gas  mixture  molecular  weight,  and  ^  is  the  molecular  weight  of  species  i.  The 
Stefan-Maxwell  relations  may  also  be  used 


where  x^  is  the  mole  fraction  of  species  i  and.^j  is  the  binary  diffusion  coef¬ 
ficient  of  species  i  into  j .  Both  of  these  expressions  are  complex  in  that  the 
multicomponent  diffusion  coefficients  are  difficult  to  evaluate,  and  the  Stefan- 
Maxwell  relations  provide  only  implicit  expressions  for  the  j^.  For  the  special 
case  when  all  diffusion  coefficients  can  be  assumed  equal  and  thermal  diffusion 
can  be  ignored,  Fick's  law  results: 
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This  technique  is  not  used  in  this  analysis.  A  different  simplification  is 
used  to  work  in  terms  of  "elemental"  conservation  rauher  than  species  conserva¬ 
tion  however.  The  term"element"  (in  quotes)  is  used  to  refer  to  tnose  atoms  or 
groupings  of  atoms  which  according  to  equilibrium  relations  are  conserved.  Ref¬ 
erence  7  discusses  the  merits  of  this  approach  in  more  detail.  Defining  as 
the  mass  fraction  of  "element"  k  in  species  i,  multiplying  the  species  equations 
by  ot^,  and  summing  over  all  species  results  in  the  following  conservation  of 
"elements"  equations: 
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where  Kk  is  the  mass  fraction  of  "element'  «  . 


it  f  .  *.e  J  by 
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It  has  also  been  assumed  that  all  =  eD.  Tne  "elemental"  approach  results  in 

significantly  fewer  simultaneous  equations  than  the  conservation  of  species  ap¬ 
proach,  and  the  equating  of  all  eD^  gives  sufficiently  accurate  solutions  for 
most  types  of  problems. 

The  streamwise  momentum  equation  can  be  written  as 
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where  P  is  the  local  static  pressure,  and  eddy  viscosity  eM  is  defined  in  terms 
of  the  Reynolds'  stresses  of  turbulent  flow  by 
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The  transverse  direction  momentum  equation  reduces  to  zero  when  longitudinal 
curvature  effects  are  ignored. 
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The  energy  equation  for  this  general  chemistry  boundary  layer  is 
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where  H  is  the  total  enthalpy  (static  plus  kinetic) 


ht  =  h  +  T 


(14) 


h  is  the  static  enthalpy  including  chemical  as  well  as  sensible  contributions 


h  = 


(15) 


h^  is  the  static  enthalpy  of  species  i 


hi 
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(16) 


T  is  the  temperature,  h?  is  the  heat  of  formation  of  spcies  i  at  the  reference 
temperature  T°,  Cp.  is  the  specific  heat  of  spcies  i,  is  the  frozen  specific 
heat  of  the  gaseous  mixture 


C 


P 


(17) 


X  is  the  thermal  conductivity,  R  is  the  gas  constant,  x^  is  the  mole  fraction  of 
t  ecies  j,  the  turbulent  enthalpy  transport  coefficient  is  defined  by 
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and  gr  is  the  net  one-dimensional  energy  flux  towards  the  surface  due  to  radia¬ 
tion  absorption  and  emission. 

In  the  energy  equation,  as  in  the  species  conservation  equations,  it  is  nec¬ 
essary  to  evaluate  molecular  diffusion  flux  j^.  As  discussed  earlier,  the  gen¬ 
eral  expressions  for  these  terms  are  difficult  to  work  with,  therefore  an  approx¬ 
imate  technique  for  multicomponent  diffusion  has  been  derived  in  reference  8  and 
is  used  in  the  present  analysis.  Since  the  present  emphasis  is  on  turbulent 
flow  problems  rather  than  molecular  diffusion,  only  the  results  of  the  approxi¬ 
mations  for  diffusion  are  presented  here.  Approximating  ^ ^  by 


and  defining 
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the  species  and  "elemental"  laminar  flux  relations  can  be  expressed  as 
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The  "elemental"  species  conservation  equation  becomes 
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while  the  energy  equation  can  be  expressed  as 
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If  equal  diffusion  coefficients  are  assumed,  p,  =  \/%,  C  =  C  ,  and  h  =  h.  When 
thermal  diffusion  is  to  be  neglected,  cfc  =  0  and  p^  =  In  p2. 

Equations  (4),  (11),  (33),  and  (34)  comprise  the  boundary  layer  conservation 
equations,  including  the  approximations  for  unequal  thermal  and  multicomponent 
diffusion  coefficients  of  reference  8.  The  equations  are  parabolic  in  nature, 
therefore  requiring  specifications  of  the  dependent  variables,  their  derivatives. 
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or  a  linear  combination  thereof  along  the  wall  (y  =  0) ,  edge  of  the  boundary 
layer,  and  at  the  initial  body  station.  Typical  sets  of  boundary  conditions 
will  be  discussed  later  in  this  report.  Also  necessary  in  the  mathematical  for¬ 
mulation  of  the  problem  is  the  specification  of  the  molecular  transport  proper¬ 
ties,  equation  of  state  and  equilibrium  (or  nonequilibrium)  relations  for  the 
multicomponent  gas,  and  a  description  of  the  eddy  viscosity,  conductivity  and 
diffusivity.  The  molecular  transport  properties,  equation  of  state,  and  equilib¬ 
rium  relations  are  discussed  in  references  2  and  7.  The  turbulent  flow  model  is 
discussed  in  the  next  subsection. 

2.  TURBULENT  FLOW  CONSIDERATIONS 

In  the  conservation  equations  developed  above,  .he  concepts  of  eddy  viscos¬ 
ity,  eddy  diffusivity,  and  eddy  conductivity  were  used  to  express  the  correla¬ 
tions  of  fluctuating  velocity,  species,  and  enthalpy  fields  in  terms  of  mean 
field  quantities.  This  is  only  one  of  several  possible  techniques  of  closing 
the  set  of  equations  (assuming  satisfactory  expressions  for  the  eddy  parameters 
are  available) ,  and  it  does  not  provide  any  information  regarding  the  evolution 
of  the  turbulent  correlations  as  the  flow  progresses  downstream.  Admittedly,  it 
would  be  more  desirable  to  describe  the  turbulent  fluctuations  in  a  more  complete 
manner  such  as  with  an  entrainment  relation,  turbulent  kinetic  energy  relation, 
or  a  local  turbulent  constitutive  equation  (reference  9) .  However,  these  tech¬ 
niques  are  still  in  early  stages  of  development  even  for  incompressible  single 
component  flows,  therefore  a  more  proven  approach  was  selected  for  the  present 
analysis.  The  Boussinesq  description  of  turbulent  boundary  layers  has  proved  to 
be  very  useful,  particularly  for  complex  reacting  flows  such  as  are  being  de¬ 
scribed  here,  and  will  be  used  exclusively  in  the  present  analysis. 

There  is  a  wide  amount  of  latitude  possible  even  within  the  eddy  viscosity 
framework  of  turbulence,  particularly  in  applying  classical  incompressible  models 
to  compressible  flows.  The  following  two  subsections  describe  how  the  turbulence 
model  described  in  references  3  and  10  was  applied  to  the  present  compressible 
flow  problem. 

a.  Wall  Region 

Following  the  work  of  Clauser  (reference  11)  the  boundary  layer  is  di¬ 
vided  into  a  law  of  the  wall  region  and  a  wake  region.  The  relatively  thin  wall 
region  of  the  turbulent  boundary  layer  is  characterized  by  very  steep  gradients 
in  the  turbulent  transport  and  mean  field  properties.  Turbulent  stress  varies 
from  zero  at  the  wall  to  near  its  maximum  value  at  the  outer  edge  of  the  wall 
region.  There  is  a  vast  amount  of  empirical  evidence  that  these  turbulent 
stresses  and  also  the  mean  flow  field  properties  can  be  described  entirely  in 
terms  of  the  wall  state,  wall  fluxes,  thermodynamic  and  transport  properties  of 
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the  fluid,  and  the  normal  coordinate  y.  Since  the  streamwise  coordinate  does  not 
enter  the  solution  for  this  region,  the  problem  becomes  a  one-dimensional  initial 
value  problem.  Eliminating  s  derivatives  from  the  continuity  equation  and  neg¬ 
lecting  variations  in  r  due  to  the  thinness  of  the  layer  results  in 


P&l  =  o 

dy 


(35) 


or 


pv  =  pwvw  (36) 

where  the  subscript  w  refers  to  the  wall  value.  Thus  the  wall  injection  rate, 
Pwvw,  which  may  be  a  function  of  s,  determines  the  transverse  mass  flux  through 
the  entire  wall  region.  Using  the  same  technique  for  the  momentum  equation  and 
substituting  equation  (36) 

pwvwU  =  P(V  +  CM>  37  "  Tw  (37) 

where  the  wall  shear,  is  also  typically  a  function  of  s.  For  flows  over  an 
impermeable  wall  with  constant  properties,  this  equation  reduces  to 

P(v+eM)g=Tw  (38) 


or 


t  =  tw  (39) 

indicating  that  shear  can  be  considered  constant  in  the  wall  region.  For  flows 
with  injection  or  ablation,  it  is  seen  that  shear  varies  with  the  mass  injection 
rate  and  local  velocity,  that  is, 

T  =  Tw  +  pwvwu  (40) 

This  one-dimensional  description  of  turbulence  in  the  wall  region  will  be  useful 
in  formulating  a  mixing  length  model  for  eddy  viscosity  as  described  in  the  fol¬ 
lowing  paragraphs.  It  should  be  made  clear  however,  that  only  the  wall  region 
turbulent  shear  stress  is  assumed  to  behave  in  a  one-dimensional  fashion.  In 
the  solution  procedure,  the  complete  two-dimensional  equations  of  motion  are 
solved  over  the  entire  boundary  layer. 
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A  complete  investigation  of  the  validity  of  the  mixing  length  postulate  for 
fJows  with  injection  has  been  reported  in  reference  12.  The  analysis  used  in 
this  investigation  is  an  extension  of  that  work;  therefore,  the  reader  should 
refer  to  reference  12  for  more  details. 

Because  of  the  current  lack  of  understanding  of  turbulent  mechanisms,  "theo¬ 
retical"  predictions  of  the  variation  of  turbulence  near  the  wall  must  rely  on 
empirical  input  into  relations  based  on  some  phenomenological  dependence.  The 
generality  of  the  ultimate  goals  of  this  analysis  and  the  desire  to  approximate 
the  physical  situation  dictated  certain  prerequisites  for  the  turbulent  transport 
relations.  These  were: 

a)  The  relations  must  indicate  a  continuous  variation  of  the  turbulent 
transport  properties  from  the  wall  to  the  fully  turbulent  region 

b)  The  relations  must  be  generally  applicable  to  mass,  momentum,  and  energy 
transport 

c)  The  relations  must  be  applicable  to  compressible  or  incompressible  flows 
with  real  gas  properties 

d)  The  relations  should  be  suitable  for  transpired  and  untranspired  bound¬ 
ary  layers  without  any,  or  a  minimum,  modification  of  form. 

Two  basic  variations  of  the  eddy  viscosity  hypothesis  have  been  proposed  in 
the  past.  The  first  type  predicts  the  variation  of  turbulent  viscosity  from  the 
wall  to  the  fully  turbulent  region.  Reichardt,  Rannie,  and  Deissler,  in  refer¬ 
ences  13  through  15,  have  proposed  such  variations.  The  second  type  of  hypothe¬ 
sis  involves  a  variation  of  mixing  length  from  the  wall  into  the  fully  turbulent 
portion  of  the  boundary  layer.  Rotta,  von  KSrmSn,  and  van  Driest  (references  16 
to  18) ,  have  adopted  this  procedure.  Data  indicate  that  surface  mass  addition 
strongly  affects  the  eddy  viscosity  profile,  and  it  was  found  that  the  first 
type  of  hypothesis  could  not  be  simply  modified  to  predict  this  variation.  On 
the  other  hand,  success  of  the  mixing  length  theory  in  predicting  profiles  in 
the  fully  turbulent  portion  of  the  boundary  layer  with  surface  mass  addition  has 
been  noted,  for  example,  in  references  19  and  20.  It  has  generally  been  concluded 
that  the  slope  of  the  linear  relation  between  mixing  length  and  distance  from  the 
wail  is  insensitive  to  surface  mass  addition.  As  a  consequence  of  this  apparent 
generality  of  the  mixing  length  approach,  it  was  adopted  for  the  present  studies. 

The  basic  mixing  length  postulate  can  be  expressed  as 

(pV>'U'  =  P*2  (§)2  =  P£M  37  (41) 

where  the  mixing  length,  1,  is  a  combination  of  various  correlations,  but  re¬ 
tains  some  relationship  to  the  scale  of  turbulence.  Prandtl  proposed  that  this 
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length  will,  in  its  simplest  form,  be  related  to  the  distance  from  a  wall,  at 
least  in  the  region  of  development  of  turbulence.  His  proposition  that 


d£ 

a?  = 


constant, 


K 


(42) 


has  been  tested  under  a  variety  of  conditions  and  found  to  be  quite  adequate  in 
the  fully  turbulent  portion  of  the  wall  region. 

As  the  wall  is  approached  however,  this  simple  relation  is  no  longer  appro¬ 
priate,  and,  in  fact,  it  can  be  shown  theoretically  that 


£im  £ 
y+o 


£im 

y-*-o 


d£ 

dy 


(43) 


This  is  a  consequence  of  the  Reichardt-Elrod  criterion  (see  reference  12) .  Thus, 
two  criteria  are  specified,  namely,  Prandtl's  hypothesis  which  is  appropriate  in 
the  fully  turbulent  portion  of  the  wall  region  and  the  Reichardt-Elrod  wall  cri¬ 
terion  as  expressed  by  equation  (43) . 

Several  means  of  expressing  a  relation  covering  the  full  range  of  y  and  in¬ 
cluding  these  limiting  criteria  have  been  used  by  other  investigators.  It  is 
advantageous  in  considering  extensions  of  mixing  length  theory  to  establish  some 
physical  logic  for  the  selected  relation.  Unfortunately,  the  understanding  of 
transition  from  the  laminar  to  the  turbulent  portions  of  the  layer  has  not 
reached  a  state  permitting  any  quantitative  specification.  Therefore,  the  se¬ 
lected  model  can  be  based  only  on  qualitative  understanding  of  the  process,  di¬ 
mensional  considerations,  and  the  above  limiting  criteria.  These  criteria  are 
satisfied  for  incompressible  flows  by  a  simple  implicit  relation  of  the  form 

||  «  (Ky  -  £)  (44) 

which  implies  that  the  rate  of  increase  of  the  mixing  length  is  proportional  to 
the  difference  between  the  value  postulated  by  Prandtl  (Ky)  and  its  actual  value. 
This  rate  of  increase  is  assumed  to  be  augmented  by  the  local  shear  and  retarded 
by  the  local  viscosity.  Using  these  parameters  to  nondimensionalize  the  above 
relation  yields 


d£ 

dy 


(Ky  -  t) 


(45) 


-14- 


AFWL-TR-69-106 


where  y*  is  the  constant  of  proportionality.  The  coefficients  K  and  y*  were 
shown  in  reference  12  to  be  invariant  for  a  wide  variety  of  flow  conditions  at 
values  of  0.44  and  11.83,  respectively. 

For  compressible  flows,  the  physical  arguments  must  be  changed  somewhat. 
Rather  than  describing  the  scale  of  a  turbulent  eddy,  it  seems  appropriate  to 
describe  the  mass  of  the  eddy,  pi,  with  respect  to  the  mass  available,  /  P  <3y. 
Thus,  by  analogy  to  equation  (44),  the  rate  of  increase  of  the  mass  of  an  eddy 
will  be  taken  to  be  proportional  to  the  difference  between  the  mass  available 
between  the  wall  and  the  point  of  interest  (times  an  appropriate  constant)  and 
the  mass  of  the  eddy: 


pdy  -  pi 


Nondimensionalizing  as  above, 


-  (k  f*  pd,  -  p.) 


/?7p 

+ 


(46) 


(47) 


The  constants  K  and  yi  are  left  at  their  incompressible  values  of  0.44  and  11.83 
for  the  time  being.  The  integral-differential  character  of  this  mixing  length 
equation  indicates  a  difficult  solution  procedure  in  the  physical  coordinate 
plane.  However,  in  the  (n,5)  coordinates  introduced  by  the  Levy-Lees  transfor¬ 
mation,  the  mixing  length  equation  simplifies  somewhat.  This  will  be  discussed 
further  in  Section  11.3. 


For  the  special  case  of  constant  properties  and  zero  injection  (constant 
shear) ,  equation  (47)  can  be  integrated  to  yield 


where 


(48) 


(49) 
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It  can  be  seen  that  the  Reichardt-Elrod  criteria  is  satisfied  at  the  wall.  For 
large  y,  Rotta's  (reference  16)  expression 


4  -  5r  <*+  -  # 

T 


is  obtained.  This  special  case  result  for  constant  property  zero  injection  flows 
is  not  used  in  the  general  analysis  technique  presented  here. 

b .  Wake  Region 

The  wake  region  of  a  turbulent  boundary  layer  is  so  named  because  the 
flow  in  this  region  tends  to  have  a  wake-like  character.  In  particular,  the 
outer  80  to  90  percent  of  the  boundary  layer  combined  with  the  local  turbulent 
eddies  dominates  the  mixing  processes  within  the  flow,  and  the  viscous  effects 
become  second  order.  Gradients  in  the  wake  region  are  typically  much  smaller 
than  those  of  the  wall  region.  Since  the  pressure  gradient  and  streamwise  deriv¬ 
ative  terms  are  important  in  the  wake  region,  the  two-dimensional  character  of 
the  turbulence  must  be  considered  in  its  entirety,  as  opposed  to  the  approxima¬ 
tions  of  the  wall  region. 

A  fortunate  feature  of  the  wake  portion  of  the  boundary  layer  is  that  eddy 
viscosity  is  nearly  constant  across  this  region,  at  least  for  equilibrium*  in¬ 
compressible  flows.  In  particular,  Clauser  (reference  11)  was  able  to  relate  the 
eddy  viscosity  to  edge  velocity  and  a  length  scale  6* 


eM  =  0.018  u^6* 


for  a  great  quantity  of  experimental  data  taken  in  equilibrium  flows. 
The  quantity  6*  in  this  relation  is  the  displacement  thickness 


in  which  the  densities  cancel  out  for  incompressible  flows.  For  compressible 
flows,  this  length  scale  is  inappropriate  since  under  some  conditions  6*  can  be 
negative.  Defining  a  velocity  displacement  thickness  as 


Equilibrium  as  used  here  refers  to  a  particular  pressure  gradient,  (5*/t  ) 
(dP/dx) ,  which  results  in  self-similar  velocity  profiles  (reference  11) . 
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“  o 

the  eddy  viscosity  in  the  wake  portion  of  the  flow  will  be  taken  as 

eM  =  0.018  u16i*  (54) 


A  satisfactory  technique  for  choosing  the  correct  eM  expression  at  any  particu¬ 
lar  body  station  is  to  use  the  wall  region  expression 


e 


M 


*2 


du 

ay 


(55) 


until  eM  exceeds  the  wake  value,  equation  (54),  at  which  point  eM  is  held  con¬ 
stant  at  the  wake  value  for  the  remainder  of  the  boundary  layer  thickness.' 

c.  Boundary  Layer  Transition 

As  can  be  seen  from  the  form  of  the  conservation  equation,  both  the 
molecular  and  turbulent  transport  terms  are  considered  simultaneously.  This  is 
necessary  since  an  accurate  description  of  the  turbulent  boundary  layer  requires 
that  the  time-averaged  fluctuation  terms  disappear  near  the  wall.  Another  rea¬ 
son  for  the  inclusion  of  these  terms  is  the  description  of  laminar  or  transi¬ 
tional  flows.  From  the  form  of  equation  (54),  it  can  be  seen  that  for  very  small 
6*  the  turbulent  stresses  will  be  small  compared  to  the  laminar  ones.  Without 
any  constraints  on  the  equations  as  stated  above,  kinematic  and  eddy  viscosities 
are  equal  at  a  velocity  displacement  thickness  Reynolds  number  of  56: 

-  _  fM  _  °’018  Ul5i 

1  V  V 


•••  “  56 

l 

This  "natural"  transition  Reynolds  number  is  too  low  for  most  situations,  there¬ 
fore  eM  is  artificially  set  to  zero  until  some  other  criterion  is  satisfied.  A 
Reynolds  number  on  momentum  thickness,  Re^,  is  currently  used  to  trigger  transi¬ 
tion.  Once  the  prespecified  transition  value  for  Refl  is  exceeded,  turbulent 
transport  properties  are  immediately  brought  into  the  solution.  Being  a  non¬ 
similar  solution,  the  influence  of  the  upstream  laminar  profile  is  felt  for  some 
distance  downstream,  thus  simulating  a  transitional  region  which  is  not  too  un¬ 
like  the  physical  situation. 
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3 .  COORDINATE  TRANSFORMATIONS 

The  equations  of  motion  for  a  boundary  layer  flow  can  be  solved  in  the  phys¬ 
ical  (s,y)  plane  by  numerous  techniques,  however  it  is  generally  advantageous  to 
transform  the  problem  to  another  coordinate  system.  The  transformed  coordinates 
offer  the  advantages  of  nondimensionalizing  the  solution,  confining  the  solution 
to  a  narrower  region,  minimizing  changes  in  the  dependent  variables,  simplifying 
boundary  conditions  and  occasionally  result  in  the  deletion  of  streamwise  deriv¬ 
ative  terms.  This  latter  possibility  occurs  only  under  very  restrictive  sets  of 
boundary  conditions.  The  coordinate  transformation  in  the  present  analysis  is  a 
variation  of  the  Levy-Lees  transformation  and  is  derived  in  its  entirety  in  Ap¬ 
pendix  I.  The  standard  Levy-Lees  transformation  takes  the  form 


/s  2  K 

plululro  ds 

o 


n 


(56) 


The  first  alteration  of  this  transformation  is  actually  a  mathematical  conven¬ 
ience  for  carrying  out  the  numerical  solution.  Introducing  a  stretching  param¬ 
eter  aH  in  the  normal  coordinate,  a  new  coordinate  system  is  defined  by 

X  -  K 

n  «  (57) 


The  parameter  aH  is  taken  as  a  function  of  X  only  and  is  determined  implicitly 
during  the  solution.  Its  purpose  is  to  stretch  n  coordinate  such  that  the 
boundary  layer  remains  of  constant  thickness  in  the  n  coordinates. 

Since  a  new  variable  c»H(5)  is  introduced,  an  additional  relation  is  required. 
This  is  conveniently  supplied  by  constraining  some  arbitrary  point  near  the  bound¬ 
ary-layer  edge,  nc»  to  have  a  specified  streamwise  velocity,  c,  near  (but  some¬ 
thing  less  than)  the  edge  value: 


f ' 


'"'edge 


(58) 
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where  f  is  the  transformed  stream  function  defined  as 


and  the  prime  denotes  differentiation  with  respect  to  n  so  that 

<60> 

Examples  of  the  utility  of  the  stretching  parameter  ccH  are  contained  in  refer¬ 
ence  2. 

The  second  change  in  the  Levy-Lees  transformation  has  to  do  with  the  trans¬ 
verse  curvature  effect.  For  very  thin  axisymmetric  bodies,  it  is  possible  to 
have  boundary  layer  thicknesses  on  the  order  of  the  body  radius  rQ.  In  this  in¬ 
stance,  it  is  necessary  to  treat  r  as  a  function  of  y,  thereby  including  its 
variation  through  the  boundary  layer.  The  coordinate  transformations  become 


5 


Pl'Wo 


ds 


(61) 


Utilization  of  the  above  coordinate  transformation  relations  results  in  a 
new  set  of  governing  equations  in  the  (|,ri)  coordinate  plane  which  will  be  given 
below.  The  hat  (A)  notation  will  be  dropped  for  the  remainder  of  the  text  for 
simplicity,  however  £  and  n  are  given  by  equation  (61).  Primes  will  refer  to 
derivatives  with  respect  to  n  except  when  noted  otherwise. 

The  global  continuity  equation  is  automatically  satisfied  by  the  definition 
of  a  transformed  stream  function  f(£,n),  shown  in  equation  (59),  and  redefined 
here  in  the  final  coordinate  system: 
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The  other  governing  equations  will  be  discussed  separately. 
Streamwise  momentum  equation 


(63) 


ff"  + 


+  6 


fi 


°_1 

P 


3f ' 

3  In  £ 


3  In  a„ 

f»2  _ £ 

3  In  £ 


f" 


if _ \ 

3  in  l) 


(64) 


In  this  equation,  utilizing  the  technique  of  reference  21,  the  transverse  curva¬ 
ture  effect  is  included  entirely  in  the  coordinate  transformation  and  in  the  def¬ 
inition  of  t: 


t  = 


1  + 


2aH/2f  cos0 


u,r 


ro 


(65) 


where  9  is  the  angle  between  the  surface  normal  and  a  plane  normal  to  the  body 
centerline  (see  Figure  1) .  Other  definitions  of  interest  are: 


C  =  -£H_ 

P^l 


(66) 


6 


3  In  u. 

2  mrr 


(67) 


For  solutions  without  consideration  of  transverse  curvature,  t  is  set  to  1.0 
throughout  the  boundary  layer. 

Turbulent  model  equations 

The  turbulent  fluctuations  are  related  to  the  mean  field  through  the  eddy 
models  described  in  equations  (2) ,  (12) ,  and  (18)  .  Eddy  viscosity  is  described 
by  a  wall  law  and  wake  law,  while  eddy  diffusivity  and  conductivity  are  related 
to  eddy  viscosity  by  turbulent  Schmidt  and  Prandtl  numbers: 
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Defining 


Sc.  =  I* 
*  eD 


eH 


6  £ 


m 

plulro 


5  =  P£ 


Re, 


PjUjfi 


p2e, 


M 


'M  -  P^ 


(68) 

(69) 

(70) 

(71) 

(72) 

(73) 


The  wall  region  eddy  viscosity  relation  becomes 


where 


eM  " 


p(Re6}  7a 


pl«H 


izf" 


(wall  region) 


0.018 


(wake  region) 


(74) 


(75) 


(76) 


Transverse  curvature  is  not  considered  in  determining  the  wake  region  length 
scale  6*.  The  governing  equation  for  mixing  length,  which  must  be  solved  for 
the  entire  boundary  layer  although  it  is  used  only  in  the  wall  region,  is 


d  i 

an 


aHp16^x7p 


V 


(77) 
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Diffusive  fluxes 


The  normalized  diffusive  energy  flux  is  given  by 


q*  =  -  — 
a  oH 


^  ui  +  i  T*  +  A  L  -  (c  +  fsLV 

1  Pr  SE  1  \  P  H»2) 


+  cfcRTp^  +  (h  -  h  +  ctRTu3)yjJ 


'K 


-|_ui  +  _E^T.  +_i_  {h,  -Cjr*> 


H 


Sc± 


(85) 


where  Pr  is  the  Prandtl  number  based  on  the  frozen  specific  heat 


Pr 


"  X 


(86) 


The  turbulent  contribution  to  the  diffusive  energy  flux  is  contained  in  the  last 
bracketed  term,  which  is  left  uncombined  with  the  other  terms  for  clarity.  The 
fact  that  the  gross  simplifications  of  the  turbulent  model  are  included  in  the 
same  equation  with  the  rather  sophisticated  unequal  molecular  diffusion  model  is 
merely  a  mathematical  convenience  stimulated  by  the  requirement  for  calculations 
in  all  types  of  flow  situations,  including  both  laminar  and  turbulent  flows.  Un¬ 
equal  molecular  diffusion  and  thermal  diffusion  effects  may  be  important  in  the 
laminar  sublayer  region  of  a  turbulent  boundary  layer,  however. 

Normalized  molecular  diffusive  flux  of  species  i  is 


_C _ 

“h5* 


(Z; 


*i»*4 


where  Sc  is  a  system  property  defined  by 


(87) 


Sc  =  -  (88) 

pDy  2 

The  Sc  is  a  Schmidt  number  based  on  the  self-diffusion  coefficient  for  a  ficti¬ 
tious  species  representative  of  the  system  as  a  whole.  The  normalized  molecular 
diffusive  flux  of  the  "elemental"  species  is 
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f5k +  (ik  -  ] 


When  certain  groupings  of  parameters  are  constant  so  that  the  flow  simi¬ 
larity  assumption  is  valid,  the  terms  on  the  right-hand  side  of  the  conservation 
equations  (equations  (64),  (79)  and  (83))  vanish,  ir.  which  case  the  conservation 
equations  become  ordinary  differential  equations.  It  should  be  emphasized  that 
the  equations  as  presented  herein  are  equivalent  to  the  corresponding  boundary- 
layer  equations  presented  in  Section  II. 1.  That  is,  no  similarity  assumptions 
have  been  made  in  their  development. 

Equations  (82),  (67),  and  (63)  for  a*,  8,  and  fw,  respectively,  are  inde¬ 
terminant  at  the  stagnation  point  of  a  blunt  body.  Special  forms  for  these  equa¬ 
tions  valid  at  the  stagnation  point  are  shown  in  reference  2  to  be  given  by 


(pwV“*  ), 


where  for  Newtonian  flow 


6SP  “  1/(K  +  1) 


=  (2P/p), 


with  R  an  effective  nose  radius  taking  into  account  the  shock  shape.  Alter- 
ef  r 

natively,  8  and  (du./ds)  can  be  computed  from  curve  fits  of  the  inviscid 
sp  -L  sp 

pressure  distribution.  The  transverse  curvature  parameter  t  also  requires  some 
special  treatment  at  a  stagnation  point.  The  troublesome  term  is  cos  6/rc 
which  is  evaluated  at  a  stagnation  point  by 

r  -»/* 


I  cos  9  \ 

l  ro  )'• 


*1  -  -s-> 
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In  addition,  to  improve  the  accuracy  of  numerical  integration  procedures 
in  the  nose  region,  %  and  f  can  be  computed  by  the  following  relations 


(95) 


-Vz 
(25)  /2 

TTTT) 


Vw  ft) 


d(sK+1) 


which  rake  advantage  of  the  fact  that  u^/s  and  rQ/s  vary  more  nearly  linearly  in 
the  stagnation  region  than  do  u^  and  r  .  These  equations  are  also  discussed 
more  thoroughly  in  reference  2. 


4 .  BOUNDARY  CONDITIONS 

The  usual  set  of  boundary  conditions  for  the  boundary  layer  flow  problem 
consists  of  the  specification  of  initial  profiles  for  the  dependent  variables 
f',  Ht,  and  K^,  plus  additional  specifications  of  these  quantities  along  the 
wall  and  at  the  edge  of  the  boundary  layer,  and  the  specification  of  fw  along 
the  wall.  However,  since  the  main  utilization  for  the  analytical  technique  pre¬ 
sented  here  is  to  compute  boundary  layer  properties  for  flows  over  ablating  or 
transpired  surfaces  (heat  shields,  rocket  nozzles,  etc.),  these  boundary  condi¬ 
tions  have  been  greatly  generalized.  The  numerous  options  resulting  from  this 
generalization  are  discussed  below. 

The  boundary  layer  edge  conditions  typically  are  found  from  an  isentropic 
expansion  from  known  elemental  gas  composition  and  stagnation  conditions.  Thus, 
given  a  set  of  stagnation  conditions  and  a  description  of  local  static  pressure 
along  the  surface  of  interest,  the  techniques  of  reference  7  may  be  used  to  es¬ 
tablish  the  entropy  of  the  gaseous  mixture  which,  when  combined  with  the  speci¬ 
fied  pressures,  can  be  used  to  establish  the  complete  equilibrium  edge  gas  state 
at  each  body  station.  Edge  boundary  conditions  then  would  consist  of 


f  *  ss  a. 

redge  aH 

^T  =  ^T  (97) 

edge  edge  actual 

\  =  h 

edge  edge  actual 
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where  the  subscript  "edge"  refers  to  conditions  specified  at  ne<J  ,  chosen  to  be 
outside  the  boundary  layer  (see  Section  II. 3).  An  additional  constraint  at  the 
boundary  layer  edge  which  is  necessary  only  when  cubics  are  used  is  the  require¬ 
ment  of  zero  slope,  i.e., 


II 


K  - 0 

edge 


(97a) 


It  is  possible  to  specify  edge  entropy  as  well  as  pressure.  The  techniques  of 
reference  7  are  then  used  to  establish  the  complete  edge  gas  state  for  a  non- 
isentropic  expansion  around  the  body  of  interest. 

Initial  profiles  of  f',  H^,  and  are  more  difficult  to  establish  for 
the  general  problem,  therefore  calculations  are  often  started  with  reasonable 
assumed  profiles  far  upstream  of  the  region  of  interest  so  that  effects  of  erron¬ 
eous  assumptions  will  die  out.  Another  possibility  for  initially  laminar  prob¬ 
lems  is  to  assume  a  similar  solution  as  a  starting  profile.  This  assumption 
reduces  the  equations  to  ordinary  differential  equations  at  the  starting  point, 
which  may  be  solved  simultaneously  for  a  set  of  profiles  unique  to  the  assumed 
edge  and  wall  state.  The  similar  solution  is  exact  at  a  body  stagnation  point, 
therefore  this  option  is  particularly  valuable  for  blunt  body  problems. 

The  wall  boundary  conditions  allow  the  widest  selection  of  options.  The 
simplest  combination  is  the  straightforward  assignment  of  velocities,  enthalpy, 
and  elemental  concentrations  at  the  wall: 


f'  =•  0 
w 

no  slip 

f  =  f  (C ) 
w  w 

specified  Pwvw 

\  -  h»m 

specified  enthalpy 
at  the  wall 

of  gas 

Kk  -  Kk  (U 
w  w 

specified  wall  gas 
composition* 

elemental 

Wa.i  temperature  may  be  used  to  find  wall  enthalpy  in  the  above  formulation. 
Also,  wall  mass  diffusive  fluxes  of  up  to  three  individual  injectants  may  be 

*  ~ 

It  is  physically  unrealistic  in  most  cases  to  assign  when  diffusion  coeffi¬ 
cients  are  unequal  since  the  contribution  to  by  preferential  diffusion  of 
the  various  "elements"  to  the  surface  is  not  known  a  priori. 
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assigned  in  lieu  of  and  ptjvw.  With  the  values  of  the  dependent  variables 
all  directly  assignee  in  this  manner,  the  boundary  layer  problem  is  uncoupled 
from  the  surface  chemistry  interaction. 

The  inclusion  of  surface  material/boundary  layer  gas  interaction  chemistry 
in  the  boundary  layer  problem  forms  the  second  major  set  of  wall  boundary  condi¬ 
tion  options.  Using  the  surface  thermochemistry  techniques  of  reference  7,  it 
is  possible  to  specify  given  mass  fluxes  of  the  (up  to)  three  injectants  at  the 
wall  and  require  chemical  equilibrium  between  the  injectants,  the  wall  material, 

and  the  adjacent  gas  stream.  In  this  instance,  the  values  of  H™  (i.e.,  T  )  and 
~  ^  w  v 

Kkw  are  found  by  simultaneous  solution  of  the  local  surface  chemical  equilibrium 

equations,  surface  mass  balances,  and  the  no-slip  velocity  boundary  conditions. 
Alternatively,  selected  chemical  reactions  at  the  wall  can  be  kinetically  con¬ 
trolled  through  Arrhenius -type  rate  law  formulations  and  included  in  the  surface 
chemistry  description. 

In  the  use  of  this  boundary  layer  technique  in  conjunction  with  in-depth 
charring  ablation  analyses,  the  chemically  active  injectants  might  result  from 
the  pyrolysis  of  an  internally  decomposing  material,  surface  material  combustion 
or  phase  change,  and  mechanical  removal.  A  variation  of  this  type  of  wall  bound¬ 
ary  condition  is  to  specify  the  wall  temperature  or  enthalpy  and  allow  the  sur¬ 
face  chemistry  calculations  to  compute  the  necessary  p  v  and  Kj^  .  In  summary, 
the  surface  equilibrium  wall  boundary  condition  is 


f;  - 

0 

no  slip 

fw  - 

fw<5> 

specified  pyvw 

ht  = 

ht 

w 

wequil 

from  surface  equi¬ 

V  - 

\ 

librium  requirement 

w 

wequil 

The  final  wall  boundary  condition  category  involves  the  use  of  a  steady 
state  energy  balance  at  the  surface.  A  general  surface  energy  balance  can  best 
be  understood  by  examination  of  a  schematic  representation  of  the  energy  fluxes 
to  an  ablating  or  nonablating  (mc  =  0)  surface: 

qa  qr  (pv)whw  ( infinitesimally  thin 

< control  volume  at 
(surface 


AFWL-TR-69-106 


Summing  terms, 


m  h 
9  9, 


w 


mh„ 
c  c. 


w 


qr  “  (pv)whw  - 


qcond 


=  0 


(100) 


which  is  valid  in  either  a  transient  or  steady-state  situation.  In  general,  an 
in-depth  charring  ablation  solution  would  be  needed  to  provide  the  conduction 
term  <3con<3  and  the  pyrolysis  gas  rate,  mg.  Under  steady  state  conditions,  the 
internal  pyrolysis  "front"  and  the  charred  surface  are  assumed  to  be  receding  at 
the  same  rate,  therefore  requiring  that  the  energy  conducted  into  the  wall  mate¬ 
rial  must  equal  the  enthalpy  rise  of  the  wall  material  and  pyrolysis  gases.  In 
equation  form 


q _ .  -  m  (h  -  h°)  -  m  (h  -  h°)  =  0  (101) 

cond  <=  cw  c  g  gw  g 

Substituting  into  equation  (100) ,  the  steady  state  energy  balance  becomes 

q  ~  q  -  (PV)Mh  +  mh°  +  m  h°  =  0  (102) 

®w  *w  w  w  c  c  9  9 

In  this  equation,  q3w  is  the  wall  value  of  the  energy  flux  defined  in  equation 
(85) ,  and  is  found  in  the  course  of  the  boundary  layer  solution.  The  surface 
equilibrium  requirement  is  always  used  in  conjunction  with  the  steady  state  en¬ 
ergy  balance.  Therefore,  if  one  specifies  the  compositions  and  heats  of  forma¬ 
tion  of  the  pyrolysis  gas  and  char  materials ,  the  simultaneous  solution  of  the 
energy  equation  above  and  the  surface  chemistry  relations  mentioned  earlier  com¬ 
pletely  couples  the  boundary  layer  flow  to  the  surface  response.  The  steady 
state  assumption  is  good  even  in  transient  situations  for  large  ablation  rates 
or  small  thermal  diffusivity  of  the  ablation  material  (reference  22) .  In  sum¬ 
mary,  the  use  of  the  steady  state  energy  balance  results  in  the  following: 


w 


0 


no  slip 


s .  s . 


steady  state 
energy  balance 


wequil 

‘kw 

wequil 


surface  equilibrium 
requirement 


(103) 
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SECTION  III 

INTEGRAL  MATRIX  SOLUTION  PROCEDURE 

The  solution  of  the  transformed  boundary  layer  equations  presented  in  Sec¬ 
tion  II  uses  an  integral  matrix  method  which  has  been  developed  specifically  for 
the  solution  of  chemically  reacting,  nonsimilar,  coupled  boundary  layers.  A  com¬ 
plete  presentation  of  the  integral  matrix  procedure  was  included  in  reference  2, 
where  solution  of  laminar  flow  problems  was  discussed.  In  the  present  effort, 
this  technique  has  remained  essentially  unchanged,  however  new  variables  and 
equations  have  been  added  to  describe  the  turbulent  aspects  of  the  flow  and  to 
include  transverse  curvature  effects.  The  present  discussion  will  therefore  re¬ 
view  only  the  highlights  of  the  method,  and  the  reader  may  refer  to  reference  2 
for  more  details. 

In  the  integral  matrix  procedure,  the  primary  dependent  variables  and 
their  derivatives  with  respect  to  n  are  related  by  Taylor  series  expansions  such 
that  these  dependent  variables  are  represented  by  connected  quadratics  or  cubics 
(either  option  is  available).  That  is,  f',  H^, ,  and  are  expanded  in  Taylor 
series  form  and  the  series  are  truncated  to  reflect  the  proper  polynomial  repre¬ 
sentation.  A  nodal  network  is  defined  through  the  boundary  layer  and  the  Taylor 
series  expansions  are  assumed  valid  between  each  set  of  nodes,  with  an  additional 
requirement  of  continuous  first  and  second  derivatives  (a  spline  fit) .  Primarily 
for  convenience,  the  conservation  equations  are  integrated  across  each  "strip" 
(between  nodal  points)  using  a  unity  weighting  function.  The  linear  Taylor  ser¬ 
ies  expansions  together  with  linear  boundary  conditions  form  a  very  sparse  ma¬ 
trix  which  has  to  be  inverted  only  once  for  a  given  problem.  The  nonlinear 
boundary  layer  equations  and  nonlinear  boundary  conditions  are  then  linearized, 
the  errors  being  driven  to  zero  using  Newton-Raphson  iteration. 

In  Section  III.l,  the  Taylor  series  expansions  are  presented,  the  inte¬ 
grated  form  of  the  momentum  equation  is  discussed,  and  techniques  for  evaluating 
integral  terms  are  demonstrated.  In  Section  III. 2,  the  special  techniques  ap¬ 
plied  to  the  mixing  length  differential  aquation  are  discussed,  and  in  Section 
III. 3  the  actual  simultaneous  equation  solution  procedure  is  summarized. 

1.  INTEGRAL  STRIP  EQUATIONS  WITH  SPLINED  INTERPOLATION  FUNCTIONS 

Consider  the  boundary  layer  in  the  region  of  a  given  streamwise  station  s 
as  being  divided  into  N-l  strips  connecting  N  nodal  points.  These  nodal  points 
are  designated  by  where  i  =  1  at  the  wall  and  N  at  the  edge  of  the  boundary 
layer.  Consider  a  function  p(n)  which  with  all  its  derivatives  is  continuous 
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in  the  neighborhood  of  the  point  n  =  Then,  for  any  value  of  a  in  this  neigh¬ 

borhood,  p(n)  may  be  expressed  in  a  Taylor  ssir ies  expansion  as 


i+l 


=  Pi  +  p!5n 


+  p'!  +  p*!'  +  p’!" 

ri  2  *16  *i 


(6n) 

24 


(104) 


where 


6n  =  ni+l  ”  ni  (105) 

Conventional  finite  difference  schemes,  in  effect,  typically  truncate  the  Taylor 
series  after  the  first  term  and  use  the  resulting  expression  to  relate  p'  to  p, 
etc . ,  that  is 


P 


i 


Pi+1  ~  pi 


<5n 


(106) 


Round-off  error  is  then  of  order  ( 6 n ) 2  and  many  nodes  must  be  chosen  to  bring 
this  value  down  to  acceptable  limits .  One  can  achieve  a  reduction  in  the  number 
of  nodes  for  a  given  accuracy  by  employing  a  quadratic  or  cubic  relation  repre¬ 
senting  the  function  p  over  the  interval  of  interest.  This  can  be  achieved  by 
truncating  the  Taylor  series  after  the  third  or  fourth  term.  The  cubic  approxi¬ 
mation  will  be  used  for  the  remainder  of  this  discussion.  The  pi  can  be  consid¬ 
ered  to  be  any  of  fi#  f|,  fV,  f|"  ,  ,  H^,.,  ,  Kk.,K£.,or  K£..  Since  the  high¬ 

est  derivatives  of  the  dependent  variables  which  appear  in  the  boundary  layer 
equations  are  f  V*  ,  &£.  and  k£.  ,  it  is  reasonable  to  truncate  the  series  at  the 
next  highest  derivative  and  to  consider  that  derivative  as  being  constant  be¬ 
tween  and  that  is, 


fllll 

i  iH 


fin  _  fin 

i+l  i 
6n 


■HI" 

1  Ti+1 


-  H  " 
_i+l  Ti 
5t) 


i+l 


K'  -  KJ1 
i+l  _i 

6t) 


(107) 


Thus,  rather  than  using  finite  difference  approximations  similar  to  equation 
(106)  which  are  substitued  directly  into  the  governing  differential  equations, 
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a  set  of  linear  relations  between  the  dependent  variables  amd  their  derivatives 
is  obtained  and  is  solved  s imu 1 taneou s ly  with  the  governing  differential  equa¬ 
tions.  These  linear  relations  are  of  the  form 


£ui -  fi +  q5” +  +  £r  +  fiVi  tt-1  - ° 


-  Pi+1  *  Pi  +  Pi«n  +  pv 


(<Sr>)  2 


i  ,  i  ,  „  5  ,  ji  6  n 

-  Pi+1  +  Pi  +  ?i  ~2  +  pi+l  T  “ 


(108) 


(109) 


(110) 


where  in  equations  (109)  and  (110)  the  represents  f|,  represents  K.x  ,  and 
represents  each  of  the  K  sets  of  Ky^ . 

Notice  that  f'  has  been  taken  to  be  a  cubic  over  each  strip,  rather  than 
the  stream  function,  f,  since  it  was  desired  to  represent  velocity  (u  =  u^f'/aH) 
with  the  cubic.  Equations  (108)  through  (110)  above,  when  written  for  each  adja¬ 
cent  pair  of  nodes,  give  (3  +  2K)(N  -  1)  simultaneous  algebraic  equations  for 
the  N(4  +  3K)  +1  unknowns,  f  ,  f  *  .  f",  f"'  ,  a„,  HT  ,  H-f  ,  H#  Kk  ,  Ki  ,  Kjl 
at  each  streamwise  station,  where  K  is  the  number  of  elemental  species.*  The 
Taylor  series  equations  are  written  for  only  K-l  species  since  the  overall  mass 
bilance  equation  supplies  the  remaining  elemental  concentration.  Additional  re¬ 
lations  must  come  from  the  governing  differential  equations  and  the  boundary 
conditions.  It  is  important  to  note  that  the  f,  f’,  etc.,  are  treated  as  indi¬ 
vidual  variables  related  by  algebraic  equations.  It  is  also  important  to  note 
that  the  coefficients  in  equations  (108)  through  (110)  are  functions  of  6n  only? 
therefore,  this  portion  of  the  resulting  matrix  need  be  inverted  only  once  for  a 
given  problem. 

The  conservation  equations  (64) ,  (79) ,  and  (83)  contain  streamwise  deriva¬ 
tive  or  "nonsimilar"  terms.  In  the  present  solution  technique,  two  or  three 
point  finite  difference  formulas  are  considered  sufficient  to  express  these  de¬ 
rivatives,  since  gradients  in  this  direction  are  not  severe.  As  in  reference  2 


r<4>  i 

[dUnTrJj 


=  d0<  )£  +  V  >  A-l  +  d2  (  h-2 


(111) 


where  (  ) refers  to  the  previous  streamwise  station, 


The  mixing  length  is  not  included  in  this  variables  count  since  mixing  length 
(as  well  as  in  the  wake  region)  is  treated  as  a  state  property. 
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do  " 


2“2-l 


JcA£-1 


d2  =  0 


for  two-point  difference  and 


(112). 


j  _  5  £A2-1  +  2A2-2 

°  £,A£-i  2A2- 2 


d  s  -  o 


2A2-2 


2A2-1  2-lA2-2 


V27a 


2  r-i 


22-2  2-1  2-^ 


(113) 


for  three-point  difference  where  typically 

2A2-1  ln  ^2  ln  ^2-1  =  ln(?2//^2-l)  (114) 

The  three-point  difference  relation  is  generally  used  unless  a  similar  solution 
is  desired  (in  which  case  dQ  *  d^  =  d2  =  0)  oj.  unless  the  point  in  question  is 
the  first  point  after  either  (1)  a  similar  solution  or  (2)  a  discontinuity  (e.g., 
where  the  body  changes  shape  abruptly,  or  where  mass  injection  is  suddenly  ter¬ 
minated)  . 

The  next  step  in  the  treatment  of  the  conservation  equations  is  their  in¬ 
tegration  across  the  boundary  layer  "strips”.  The  primary  reason  for  this  inte¬ 
gration  is  to  simplify  the  n-derivative  terms  in  the  energy  and  species  conser¬ 
vation  equations,  since  it  is  not  convenient  to  express  the  complex  q*  and  j* 
terms  in  derivative  form.  The  solution  can  actually  proceed  very  nicely  with¬ 
out  integrating  across  strips  (see  reference  10)  without  any  noticeable  change 
in  speed,  accuracy,  or  stability  for  simplified  problems  such  as  incompressible, 
nonreacting  f.v  -a.  The  weighting  function  for  integration  between  nodes  in  this 
integral  method  is  unity.  In  the  terminology  of  the  general  method  of  integral 
relations,  where  integrals  are  carried  from  0  to  °°  in  n  (reference  23),  a  square 
wave  weighting  function  is  used  which  is  unity  across  the  strip  in  question  and 
zero  elsewhe-.e.  The  equations  are  then  integrated  N-l  times  with  the  equate 
wave  applied  to  each  strip  in  succession.  Using  the  momentum  ea'iacion  as  an  ex¬ 
ample,  the  integration  from  i-1  to  i  results  in 
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•'i-l 


ff'dn  + 


t(C  +  eM) 


f” 


1  z*1  pi  /•* 

+  /  -y  dn  -  6  /  f ' 2dn 

i-l  •'i-l  *'i-i 


■t 


f  ’  (d0f  +  dif^_i  +  d2f£-2)dn 


l 

-f  f,2[doln 

•'i-l 


+  di(ln  aH»l-l 


+  d2f4-2)dn 


+  42<ln  ■,!.>  -  Z  f«o£  + 

•'i-l 


dlf£-l 


(115) 


The  Taylor  series  approximations  introduced  earlier  can  also  be  used  to  express 

the  integral  terms  above.  As  demonstrated  in  reference  2,  the  term  f1  f’p  dn 

'i-l 


becomes 


where 


r 

I  £’  p  dr)  -  f!  xp,  +  ft’  xp.  +  f‘."xp_  +  f\" 
J  1  1  1  2  i  3  1-1 


XP 


XP, 

X 


1-  (fi  -  Pi  ^  *  Pi  ^  *  Pi-!  -rr1) 

xp3  .  «,)*(^.  -  Pi  ♦  pr  Mf.) 

“4  -  «"'•(&  -  Pi  OT  *  Pi  +  Pi-! 


(116) 


(117) 


This  technique  is  used  to  rewrite  each  of  the  integral  terms  in  equation  (115) 
above  of  the  form  f'p  dn.  The  remaining  integral  term  in  the  momentum 

equation,  J*  *  (p^pjdn,  and  the  "elemental"  source  term  in  the  "elemental"  con¬ 
servation  equation  are  evaluated  by  approximating  these  functions  as  cubics  over 
the  atrip  and  integrating  directly.  This  yields 

L  Pi 

T  dn 

-1 
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Similarly  for  the  integral  of 


/ 


♦k  dn 


(119) 


These  approximations  are  not  quite  as  good  as  the  approximations  for  f’,  HT  and 
since  continuity  of  derivatives  is  not  guaranteed  at  the  nodal  point. 

Direct  substitution  of  these  approximations  for  integral  terms  into  the 
governing  equations  results  in  the  following  forms. 


Momentum 


H 


:■(, 


t  (C  +  £  ) 

a.." 


i-1 


+  *  <4 


,pl  pl  \«n  .  /pr''i  pl pi-l\ (6n) 

,pi  +  pi-ij  2  l“iT  ■  ~iX]^ 


1  +  P  +■  dr 


dlaHi_1  +  d2aH/„2 
ttH 


f!  XP.  +  f  **  XP. 
1  1  1  2 


Energy 


+  fi  XP3  +  fi-lXP4 


Piefi 


f jZP1  +  fVZP2  +  f V'  ZP3 


+  f'."  ,ZP. 
l-l  4 


=  0 


p .  =f ! 
l 


(120) 


t(-  q*  +  q£)  +  Ht  ((1  +  do)  f  +  d1fJl_1  +  d2f 


i-1 


-  (1  +  2do) 


fiXPl  +  fiXP2  +  fi  XP3  +  f i'-lXP4 


JPi=HT. 


fiZPl  +  fiZP2  +  f i'  ZP3  +  fi-lZP4 


Pi=HT. 


Ht  ZP1  +  H£  ZP2  +  H"_ZP3  +  H"  Zl 


=  0 


p.=f ! 

*1  l 


(121) 
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"Elemental"  Species 


E^L  K,'  -  +  K.  /(I  +  dj  f  +  d.  f.  .  +  d-  f.  0 

cHsct  ‘  o  1  2  ‘-2 

+  “H  [(\  +  \.,)f  -  (\  -  \_JtT 


-  (!  +  2ty  [fi  ”l  *  fi'  ^2  +  fi"  » 3  ♦  fi"i  xp„] 

-  [£i  ZP1  *  fi  ZP2  +  fi"  ZP3  +  *1-1  ™i]  P  .  ^ 

-  [\  ZP1  +  %  ZP2  +  X  ZP3  +  Si'..1  ZP4]  .  t. 


pi  ■Kk. 


(122) 


The  following  definitions  are  necessary: 


zp1  =  6n 


yp2  T-  +  YP3  ^T1  +  YP4 


zp2  = 


(^)2(-ii  -  yp2  t  +  YP3  iinir1  +  YP4  tt") 


zp3  =  (6n ) 3 1 


116n  .  w  11  (6n)2  ,  „„  5(6n)2\ 

YP2  -"120 ~  +  YF3  “ W  +  YP4  “W/ 


,/YP, 


!!l  .  vp  in  *  YP  +  YP  d")*' 

xn  YP2  30  Yl3  HHTJ  4  252  t 


(123) 


YP1  dlp£-l,i  +  d2p£-2,i 
Yp2  =  ^2.P£-1 ,  i  +  d2P£-2,i 
YP3  =  dlp£-l , i  +  d2P£-2 , i 
YP4  =  dip£-l  i-l  +  d2p£-2  i-1 


(124) 
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and  is  defined  adjacent  to  the  brackets  in  each  term  that  uses  these  defini¬ 
tions  . 

The  conservation  equations  provide  (K+l) (N-l)  more  equations  for  the 
N(3K  +  4)  +  1  unknowns,  thereby  closing  the  problem.  However,  before  discussing 
how  this  set  of  algebraic  equations  is  solved,  Section  III. 2  describes  in  detail 
how  the  mixing  length  differential  equation  is  solved. 

2.  SOLUTION  OF  THE  MIXING  LENGTH  EQUATION 

The  mixing  length  equation  is  a  first  order  linear  differential  equation 
whose  solution  can  be  written  directly  in  general  terms.  The  differential  equa¬ 
tion  is 


dl 

d" 


(K«Hn  -l) 


(77) 


Defining 


P(n) 


(XjjP^/tTp 


(125) 


results  in 

-  <x<V  -  Dp 

The  solution  to  this  equation  is 


The  remaining  problem  is  to  evaluate  the  integral  terms.  Defining 


(126) 


(127) 


(128) 
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yields 

l  =  KoH(n  -  L)  (129) 

Appendix  II  presents  a  complete  description  of  the  technique  used  to  evaluate 
L(n) .  In  essence,  P(n)  is  assumed  to  vary  linearly  over  the  interval  to 

r,^,  and  the  integrals  are  expressed  in  a  more  tractable  form.  The  final  expres¬ 
sion  is 


The  Dawson  Integral,  Dw(  ),  can  be  evaluated  from  tables  (reference  24)  or  by  a 
series  method.  A  series  evaluation  method  is  used  in  the  present  analysis.  Thus, 
combining  equations  (129)  and  (130) ,  an  explicit  recursion  formula  for  mixing 
length  at  each  node  is  obtained.  This  mixing  length  is  a  function  of  local  shear, 
viscosity,  and  density  through  the  variation  of  P(n),  and  is  re-evaluated  at 
each  node  on  each  iteration  during  the  course  of  a  solution. 


3.  NEWTON-RAPHSON  ITERATION  FOR  A  SOLUTION 

A  complete  description  of  the  Newton-Raphson  iteration  procedure  as  ap¬ 
plied  to  the  laminar  equations  of  motion  was  given  in  reference  2.  Since  the 
procedure  is  basically  unchanged  with  the  addition  of  turbulent  flow  and  trans¬ 
verse  curvature  equations,  it  will  be  reviewed  only  briefly  here,  with  emphasis 
on  the  recent  additions . 
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To  illustrate  the  Newton-Raphson  method,  consider  two  simultaneous  non¬ 
linear  algebraic  equations 


F(x,y)  =  0  G(x,y)  =  0 


(135) 


the  solution  for  which  is  given  by  x  =  x,  y  =  y.  Define  x„  and  y  as  the  values 
of  x  and  y  for  the  m  iteration.  The  desired  solution  f(x,y)  can  be  expressed 
in  a  Taylor  series  expansion 


-  -  _  3F(x  ) 

F(x,y)  =  F (x  ,y  )  +  (x  -  x  )  - r2 — 1 5L 

m  m  m  3x 


3F(x  ,y  ) 
+  {y  -  v  )  - IS — SL 

y  ym;  3y 


+  •  •  . 


oG  (x  ,y  ) 

G(x,y)  =  G(xm,ym)  +  (x  -  xm)  m. 


(136) 


3G(x  ,y  ) 

+  (y  -  y  )  - — — 

1  ym  3y 


+  ... 


The  Newton-Raphson  method  consists  of  replacing  (x,y)  by  (xm+1,ym+1)  on  the  right- 
hand  side  of  these  expressions  and  neglecting  nonlinear  terms  in  xm+1  -  xm  and 
ym+i  -  ym.  This  yields  the  set  of  simultaneous  equations 


3F(Vym}  .  .  3F(VV 

Axm  - 3^ -  +  Aym  — - -  "  F  xm'ym 


av  +  Av  3G(Xm>ym>  =  _  . 

xm  3x  ym  3y  ”  Xm,ym 


rm  3y 


(137) 


or  in  matrix  form 


3F(x  ,y  ) 
m  m 


3F  ( x  ,y  ) 
m  m 


-  F(x  , v  ) 
m  rm' 


3G(x  ,y  ) 
m  1  m 


3G(x  ,y  ) 
m'-Mn 


-  G (x  ,y  ) 
m  m 
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where 


Ax  =  x  . , 
m  m+1 


^m  H 


(139) 


The  Ax  and  Ay  are  the  corrections  to  be  added  to  x  and  y  ,  respectively,  to 

yield  the  values  of  the  dependent  variables  for  the  m  +  lr  iteration.  Here 

F(x  ,y  )  and  G(x  ,y  )  are  the  values  of  the  original  functions  F(x,y)  and  G(x,y) 
m  m  m  m 

evaluated  for  x  =  x  and  y  =  y  •  As  the  corrections  approach  zero,  the  F{x  ,y  ) 

m  Jm  m  in 

and  G(xm,ym)  approach  zero.  Hence,  it  is  appropriate  to  look  upon  these  as  errors 
associated  with  the  original  equation  (135).  It  is  apparent  that  this  procedure 
can  be  extended  to  an  arbitrary  number  of  functions  and  a  corresponding  number  of 
primary  variables. 

For  the  purpose  of  the  present  analysis,  it  has  been  found  most  convenient 
to  consider  the  primary  variables  as  fi,  fj,  fV,  f V'  ,  H-p^,  H^,  H^,  Kj^,  , 
K£.,  and  c*H.  This  amounts  to  (3K  +  4)N  +  1  unknowns  where  N  is  the  number  of 
nodes  and  K  is  the  number  of  elemental  species  to  be  considered  in  the  boundary 
layer.  Recounting  the  number  of  equations,  we  have 


Taylor  series  expansions 
Boundary  layer  equations 
Boundary  conditions 


aH  definition 


Eqn.  Numbers 


(103)  -  (110) 

(120)  -  (122) 

(97) , (97a) ,  (98) 
or  equivalent 


Total 


No.  of  Equations 


(N  -  1)  [5  +  2  (K  -  1)  ] 
(N  -  1)  (K  +  1) 

3K  +  4 


N(3K  +  4)  +  1 


Other  secondary  variables  such  as  e,  p,  T,  etc.  are  expressed  in  terms  of  those 
listed  above.  The  corrections  in  these  secondary  variables  are  therefore  found 
in  terms  of  the  corrections  to  the  primary  variables. 

The  use  of  the  Newton-Rapbson  technique  for  the  current  set  of  equations 
requires  the  evaluation  of  the  partial  derivatives  of  each  equation  with  respect 
to  each  variable.  The  partial  derivatives  of  the  Taylor  series  equations  and 
linear  boundary  conditions  are  exactly  the  same  as  in  reference  2.  The  deriva¬ 
tives  of  the  conservation  equations  are: 

Momentum 


t(C  +  eM)f"  fAf"  AC  AeM  AotH  At  \  n  .  ,  .  i  r  ,  i  r 

~ h —  +  -  *  •  *r  T  11  a°,f  +  4£ 


+  £'  (1  +  dol4£  -  6“h  h  1  *  S-T  -  f 
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/  o.  V 


'I* 


i  1  -  ^  ^±1  Ap;  •  +  &  AP-  L 


)  Ap-  .  +  ^2-  ip!  ,  I  }  +  8a  6n  ^  1  +  -r— 

3  p-.J  l-x  *  1-1  /  H  p.  Pi_x 


6r>i  pi  Pi-1 


a,  a„  +  d-ctu 
^Jt-l  2  "  4,-2 


1+B* 


x  f  j  AXPL  +  fV  AXP2  +  fV'  ^xp3  +  fV'^  AXP4  +  XP1Af[  +  XP2Afi 


+  XP,Af V  +  XP  .  Af  V 
3  l  4  i 


i  /diCVi +  d2<v2\  r 

*-^,.*1  '  v — [qxPi  +  fp 


- 

+  fV'  XP,  +  fV'  -XP.  Aau  -  2  ZP.Af!  +  ZP-Af'.’  +  ZP,Af'." 

1  J  x~x  *t  n  xx  xx  Jx 

Jp,=f!  L 


+  ZP  .  Af  V'  - 

H  1-1 


=  -  ERROR 


Pi=fi 


where  the  ERROR  is  given  by  the  left-hand  side  of  equation  (120)  evaluated  for 
fch 

the  m  iteration. 


t(-Aq*  +  Aq*)  +  (-q*  +  q*)At  +  ^(1  +  dQ)  f  +  djf^j.  +  d2fz~2 ^ 
+  Ht(1  +  d0)Af  ^  -  (1  +  2d0)£f!  AXPX  +  fV  AXP2  +  f V'  AX?3 


(1  +  2dQ)  f|  AXP1  +  f|  AXP2  +  f V'  AXP3 


+  fV^  AXP4  +  XP-jAfj  4  XP2AfV  +XP3AfV'+  XP4AfVL}]  -fzP-jAfj 

Pi=Hx.  L 


+  ZP0A  f  V  +  ZP,Af!."  +  z.p.Af; 

XX  J  l  ■»  J 


Pi=HT. 


-  sp1aht 


+  zp2ah; 

*i 


+  Zp3AHT.  +  ZP4AHT,  ,1 

1  1  ~  J  P  i =  f 


=  -  ERROR 


(141) 
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where  the  ERROR  is  given  by  the  left-hand  side  of  equation  (121)  for  the  mt*1 
iteration  and  Aq*  is  given  by 

3 . 
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on 

6 


6jn 

6 


AaH 


-  (1  +  2d0) 


ffAXP,  +  fVAXP,  +  f  **’  AXP,  +  fV  .  AXP .  +  XP.Af! 
i  1  i  2  l  j  l-l  4  1  l 


+  XP2AfV  +  XPjAfV*  +  XP4AfV‘_1 


J?i=Kk. 


ZP ,  A f  !  +  ZP-Af 

1  l  2  l 


+  ZF3AfV' 


ZP4Af 


r-il 


'Pi-s 


ZPj^AK,  +  ZP2AK'  +  ZP3AK" 
Ki  i  i 


zp4ak" 

4  Ki-1 


=  -  ERROR 


Pi-fi 


(143) 


where  the  ERROR  is  given  by  the  left-hand  side  of  equation  (122)  evaluated  for 
the  mth  iteration  and  Aj£  is  given  by 


+ 


-  Kk)p; 


_  ASc\ 
H  Sc  / 


+  AZk  l-  (Zk  -  Kk)Zu^  +  li^(AZk  -  AKk) 


(144) 


The  technique  of  relating  corrections  on  secondary  variables  such  as  C, 
T,  Pr,  etc.,  to  corrections  in  primary  variables  was  fully  explained  in  ref¬ 
erence  2.  The  same  techniques  are  used  for  the  new  corrections  At  and  As... 

M 

Once  the  correction  coefficients  (partial  derivatives  with  respect  to 
each  primary  variable)  for  each  equation  at  each  nodal  point  are  found,  they 
are  arranged  in  matrix  form  for  further  manipulation.  The  order  of  the  primary 
variables  and  the  order  of  the  equations  is  of  some  importance  in  the  matrix 
formulation.  It  is  most  convenient  to  divide  the  variables  into  "linear" 
(symbol  L)  and  "nonlinear"  (symbol  NL)  sets,  namely 


[al_  : 

BL  ] 

”'VL  ' 

EL 

[anl  ! 

BNLJ 

AVNL 

ENL 
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where  the  linear  equations  are  the  Taylor  series  equations  and  some  of  the  bound¬ 
ary  conditions.  The  purpose  of  the  partitioning  is  to  allow  operations  on  sec¬ 
tions  of  the  coefficient  matrix  which  result  in  significant  simplification  of 
the  overall  inversion.  In  particular,  since  the  coefficients  of  the  linear  equa¬ 
tions  are  all  constant  or  functions  of  the  fixed  nodal  spacing,  this  portion  of 
the  matrix  (the  AL  portion)  can  be  diagonalized  once  and  for  all  in  any  given 
problem.  In  essence,  the  corrections  on  the  linear  variables  AVL  are  always  ex¬ 
pressed  in  terms  of  the  nonlinear  variable  corrections  AVNL.  The  choice  of 
linear  and  nonlinear  labels  for  the  variables  is  somewhat  arbitrary,  but  care 
must  be  taken  that  the  A1  matrix  not  be  singular.  It  has  been  found  convenient 
to  arrange  the  variables  into  the  linear  and  nonlinear  groups  as  follows: 
AVLp(Af2,  Af  ^ ,  *  >  > ,  Afn,  AfJ,  Af^,...  Af£,  Af£'  ,  Af  "•  , . . .  Af£  )  ;  AVL^  (AH^, 


AHm 


*H*,, 


AH4 


AH# 


AH#  , 


AH, 


#  );  and  K-l  sets  of  AVLR  (AKj^,  AK^, 


AK 


A2  a3;  -  w  2-  n  *  n  * 

K£  , . . .  AK£  ,  AK£^,  AK£^,...  ^Kk  nonlinear  variables  are  then  arranged 


*3  n  2  n 

in  the  following  order:  AVNLp  (A«H,  Afw,  Af^,  Af#,  Af#,...  Af^) ;  AVNLy  (AH^W» 

^hTw'  ^HTn_i^  '  anc*  set3  AVNLk  (  ^kw'  ^2'*’*  ^n-l^ 

The  order  of  the  linear  equations  in  the  present  matrix  procedure  is: 


No.  of 
Equations 

3N-2 


2N 


(K  -  1)  (2N) 


Description  of  Equations 

Linear  boundary  conditions  and 
Taylor  series  for  f,  f’,  f",  f"1 

Linear  boundary  conditions  and 
Taylor  series  for  HT,  H#,  H# 

Linear  boundary  condi tions_and 
Taylor  series  for  Kk,  K£,  K£ 


The  nonlinear  equations  are  sequenced  as  follows: 


No.  of 
Equations 


Description  of  Equations 

Nonlinear  boundary  conditions 
and  a  constraint 

li 

Momentum  equation  for  each  pair 
of  nodes 


(K  -1) IN) 


Energy  equation  for  each  pair  of 
nodes  plus  wall  enthalpy  equation 

K-l  sets  of  "elemental"  species 
equations  for  each  pair  of  nodes 
plus  wall  species  equation 


Special  logic  has  been  written  for  the  matrix  inversion,  taking  advantage 
of  the  regular  sparseness  of  the  matrix.  Once  the  corrections  for  thi  linear 
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and  nonlinear  variables  are  found,  these  corrections  are  added  to  the  variaoles 
to  form  the  new  guesses.  The  magnitude  of  the  errors  for  each  equation  are 
checked  and  the  procedure  advances  to  the  next  iteration  if  the  absolute  values 
of  the  errors  exceed  prescribed  upper  limits.  If  the  errors  are  acceptable, 
iteration  is  completed  for  the  current  streamwise  position  ?.  Typically,  three 
to  six  iterations  are  required  to  reach  a  satisfactory  solution. 
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SECTION  IV 

SOME  RESULTS  FOR  MULTICOMPONENT  BOUNDARY  LAYERS 

The  fact  that  the  BLIMP  computer  program  works  well  for  a  myriad  of  multi- 
component  reacting  laminar  flow  situations  including  stagnation  points,  trans¬ 
piration  cooled  surfaces  and  rocket  nozzles  is  well  documented  in  references  25 
through  27,  among  others.  The  purpose  of  this  section  is  to  present  results  for 
the  recent  additions  of  turbulent  flow  and  transverse  curvature. 

1.  ABLATING  FLAT  PLATE  IN  TURBULENT  FLOW 

This  sample  problem  is  the  formulation  of  a  typical  turbulent  channel  or 
turbulent  pipe  flow  situation  such  as  might  be  found  in  some  of  the  major  test 
facilities  around  the  country.  The  flow  stagnation  conditions  were  PQ  =  43.4 
atmospheres,  HQ  =  2100  Btu/lb.  The  flat  plate  model  was  taken  to  be  constructed 
of  graphi'  -  phenolic,  and  the  wall  temperature  was  assigned  at  Tw  =  4760°R. 
Assuming  that  the  plate  ablates  in  a  steady  state  mode, the  chemical  composition 
of  the  virgin  material  exactly  equals  the  chemical  composition  of  the  char  plus 
pyrolysis  gas  ablation  products.  The  chemical  composition  of  graphite  phenolic 
which  was  used  is 


0.9236  lbs.  C/lb. 

0.0209  lbs.  H/lb. 

0.0554  lbs.  O/lb. 

These  numbers  were  calculated  from  data  on  graphite  phenolic  given  in  reference 
28.  Pressure  and  ablation  rate  or  mass  flux  were  also  input  to  the  BLIMP  pro¬ 
gram.  The  variations  of  these  quantities  in  the  streamwise  direction  are  shown 
in  Figures  2  and  3.  The  decision  to  specify  both  the  wall  temperature  and  abla¬ 
tion  rate  precluded  the  possibility  of  using  either  the  surface  equilibrium  or 
steady  state  energy  balance  logic  in  the  BLIMP  program.  Indeed,  no  surface  chem¬ 
istry  is  needed  at  all  when  Tw  and  m  are  specified.  This  choice  was  made  how¬ 
ever,  since  it  was  felt  that  a  more  representative  boundary  layer  would  be  ob¬ 
tained  by  foregoing  surface  equilibrium  than  by  demanding  equilibrium  and  attempt¬ 
ing  to  calculate  Tw  in  the  sensitive  diffusion  controlled  ablation  region  (see 
Section  V) . 

Using  the  above  information,  the  calculation  was  started  at  S  =  0.03125 
feet.  At  the  first  station,  streamwise  derivative  terms  are  automatically 
dropped  since  no  upstream  information  is  available.  The  solution  was  turbulent 
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over  the  entire  length  of  the  2.40  foot  ablating  plate.  Additional  solutions 
were  found  at  S  =  0.3646,  0.6979,  1.1562,  1.6979  and  2.3646  feet.  The  variations 
of  calculated  local  drag  coefficient,  Reynolds  number  on  momentum  thickness,  and 
shape  factor  may  be  seen  in  Figures  4,  5,  and  6.  The  flow  is  actually  transi¬ 
tional  in  nature,  being  weakly  turbulent  at  the  beginning  of  the  plate  and 
strongly  turbulent  at  the  end.  This  fact  combined  with  the  free  stream  accelera¬ 
tion  and  decreasing  rate  of  blowing  along  the  plate  produces  the  unusual  Re0  and 
shape  factor  variations. 

Velocity  profiles  at  three  representative  stations  may  be  seen  in  Figure 
7.  These  curves  clearly  show  the  transitional  nature  of  the  flow,  with  a  nearly 
laminar  velocity  profile  shape  near  the  beginning  of  the  plate  and  a  more  char¬ 
acteristic  turbulent  profile  near  the  end  of  the  plate.  Figure  8  presents  the 
velocity  profile  at  S  =  0.03125  feet  again,  this  time  in  semi-logarithmic  coor¬ 
dinates.  The  laminar  sublayer,  transitional  wall  region,  and  the  wake  region 
are  clearly  visible  in  these  coordinates,  with  the  data  points  shown  actually 
representing  solution  points  or  nodes  in  the  computer  solution.  This  smooth  vari¬ 
ation  of  the  velocity  ratio  from  the  wall  to  the  wake  region  with  this  solution 
technique  is  of  particular  interest  in  this  figure.  Eddy  viscosity  normalized 
by  the  edge  value  of  kinematic  viscosity  is  shown  in  Figure  9.  The  eddy  viscos¬ 
ity  is  seen  to  decrease  to  values  far  below  the  molecular  viscosity  as  the  wall 
is  approached.  Chemical  species  mole  fraction  profiles  at  three  stations  are 
shown  in  Figures  ID,  11,  and  12.  Only  the  major  species  distributions  are  shown 
in  these  figures,  although  a  total  of  42  species  were  considered.  Table  I  below 
lists  the  species  that  were  considered  and  their  maximum  concentration  in  mole 
fraction  at  the  S  =  0.6979  foot  station. 

Total  central  processor  run  time  for  this  problem  vith  13  nodes,  6  body 
stations,  and  42  chemical  species  including  4  elements  was  419  seconds  on  a  CDC 
6600. 


2.  SPHERE-CONE  CONFIGURATION  WITH  LAMINAR  AND  TURBULENT  FLOW 

A  second  sample  problem  was  selected  which  demonstrates  some  of  the  flex¬ 
ibility  of  the  BLIMP  program.  The  configuration  chosen  was  a  0.500  inch  nose 
radius,  7.5  degree  half  angle  sphere-cone  consisting  of  three  surface  materials: 
graphite,  pyrolytic  boron  nitride,  and  phenolic  carbon.  Stagnation  conditions 
were  representative  of  a  single  time  in  a  severe  reentry  trajectory  with  PQ  = 

242  atmospheres,  HQ  =  5520  Btu/lb.  Figure  13  is  a  schematic  of  the  sphere-cone 
configuration.  A  total  surface  running  length  of  5.0164  feet  was  analyzed  in 
thi3  problem,  which  required  29  body  stations  and  thirteen  nodes  through  the 
boundary  layer.  Figure  14  shows  a  portion  of  the  pressure  distribution  which 
was  assumed,  the  remainder  of  the  running  length  being  deleted  from  the 
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Figure  9.  Eddy  Viscosity  Profiles  for  Flow  Over  an  Ablating  Flat  Plate 
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Table  I 

SPECIES  CONSIDERED  IN  GRAPHITE  PHENOLIC 
FLAT  PLATE  PROBLEM 


f 


I 


Maximum  Mole  Fraction 

Species 

at  S  =  0.6979  ft. 

3 

CO 

H 

N2 

e- 

C  (s) 

°2 

C 

CH 

CHN 

CHO 

ch2 

ch3 

ch4 

CN 

CCU 


C2H 

C2H2 

c2n2 

C3H 

C3H2 

C3H3 

HN 

HO 

H2 

h2n 

h2o 


4.700 

X 

10-5* 

3.405 

X 

10"1* 

-2 

2.037 

X 

10 

7.774 

X 

ltT1 

-8 

5.084 

X 

10 

0. 

2.009 

X 

io"1 

1.487 

X 

.6* 

10 

-7* 

3.217 

X 

10 

5.546 

X 

io'2* 

6.167 

X 

10-6 
a  -  8  * 

9.211 

X 

10 

-6  * 

3.088 

y. 

10 

,  -7* 

4.455 

X 

10 

-4* 

2.667 

X 

10 
„  -2 

6.358 

X 

10 

„  -6  * 

1.898 

X 

10 

-3  * 

1.698 

X 

10 

,  -3  * 

6.622 

X 

10 

,  -3  * 

1.392 

X 

10 

_  -3  * 

3.903 

X 

10 

_  5  * 

4.818 

X 

10 

-6  * 

2.928 

X 

10 

_6 

5.014 

X 

10 

•  ^  -  2 

1.355 

X 

10 

.  -2* 

1.868 

X 

10 

2.447 

X 

io"7 

-2 

1.184 

X 

10 

53- 
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Table  I (concluded) 


Species 

Maximum  Mole  Fraction 
at  S  =  0.6979  ft. 

N 

3.511 

X 

10"  5 

NO 

3.211 

X 

io"2 

0 

2.509 

X 

io"2 

C4 

C5 

2.404 

1.427 

X 

X 

10  7 

-6* 

10  6 

c+ 

5.872 

X 

10"14* 

N+ 

1.362 
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presentation  here  in  order  to  expand  the  region  of  most  interest.  The  resulting 
8  distribution  may  be  seen  in  Figure  15.  It  is  apparent  that  a  great  deal  of 
care  must  be  exercised  in  choosing  the  pressure  distribution  data  points,  since 
8  is  calculated  from  the  pressure  gradient  implicit  in  that  data. 

Table  II  summarizes  the  vehicle  configuration  in  more  detail.  For  this 
sample  problem  the  option  of  a  quasi-steady  energy  balance  at  the  wall  was  used, 
in  addition  to  demanding  surface  equilibrium,  for  all  three  surface  materials. 
The  program  then  calculates  its  own  ablation  rate,  wall  species  concentrations, 
and  temperature  levels  based  on  the  energy  balance  and  chemical  equilibrium  re¬ 
quirements.  The  quasi-steady  ablation  mode  was  selected  to  allow  specification 
of  a  minimum  of  boundary  condition  information  and  to  demonstrate  the  use  of 
this  option  as  opposed  to  arriving  at  an  accurate  ablation  prediction.  Indeed, 
the  emphasis  in  the  current  program  development  has  been  the  incorporation  of  a 
turbulent  model  into  an  existing  laminar  boundary  layer  code.  The  selection  of 
an  ablation  model  and  the  corresponding  boundary  conditions  for  a  given  surface 
material  is  a  difficult  problem  in  itself,  requiring  careful  examination  of  the 
possible  condensed  phase  products  of  reaction  at  the  surface,  kinetically  con¬ 
trolled  surface  reactions,  mechanical  failure,  interaction  with  pyrolysis  gases, 
etc.  (reference  29).  For  these  sample  problems,  the  intent  was  to  de-emphasize 
this  procedure  in  order  to  concentrate  on  the  actual  boundary  layer  behavior. 

The  calculation  was  started  at  the  stagnation  point  and  allowed  to  transist  at 
ReQ  =  250.  Figure  16  illustrates  the  distribution  of  Cf/2  over  the  first  0.30 
feet  of  surface  running  length.  As  can  be  seen  in  the  Cf/2  distribution,  transi 
tion,  occurred  between  S  =  0.0268  and  S  =  0.0323  feet.  Reynolds  number  on  momen 
turn  thickness  and  shape  factor  streamwise  distributions  may  be  seen  in  Figures 
17  and  18.  Figure  19  illustrates  the  wall  temperature  distribution  that  results 
from  the  steady  state  energy  balance  and  surface  equilibrium  assumption.  Figure 
20  illustrates  the  calculated  quasi-steady  surface  recession  rate. 

Profile  information  for  this  problem  is  particularly  interesting  since 
both  boundary  layer  transition  and  surface  material  changes  occur  over  the  for¬ 
ward  portion  of  the  body.  Velocity  profiles  are  presented  in  Figure  21.  The 
body  stations  selected  include  the  stagnation  point  (0.0  ft.),  flow  just  ahead 
of  transition  (0.02681  ft.),  flow  just  after  transition  which  happens  to  be  over 
the  boron  nitride  just  past  the  C-BN  discontinuity  (0.03231  ft.),  and  flow  over 
the  phenolic  carbon  just  past  the  BN-phenolic  carbon  discontinuity  (0.06273  ft.) 
The  first  two  profiles  are  clearly  laminar  in  nature,  whereas  the  third  is  more 
transitional.  The  last  profile  appears  to  have  the  shape  of  a  fully  turbulent 
flow.  This  gradual  evolution  of  a  turbulent  shape  occurs  because  both  laminar 
and  turbulent  transport  terms  are  retained  in  the  equations  of  motion  (once  the 
Re0  criterion  has  been  satisfied)  and  because  the  nonsimilar  analysis  retains 
the  upstream  "history"  of  the  flow.  Species  concentration  profiles  for  the  same 
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Table  IX 

SPHERE-CONE  CONFIGURATION  SURFACE  MATERIALS  DETAILS 


Surface  Running  Length 
_ (feet)  _ 

0. 

0.004174 

0.008390 

0.012695 

0.017146 

0.021817 

0.026813 

0.032308 

0.038637 

0.046657 

0.059975 

0.062727 

0.070000 

0.075000 

0.077500 

0.087500 

0.100000 

0.125000 

0.150000 

0.200000 

0.300000 

0.500000 

0.700000 

1.000000 

1.500000 

2.000000 

3.000000 

4.000000 

5.016400 


Surface  Material 


graphite 

graphite 

graphite 

graphite 

graphite 

graphite 

graphite 

pyrolytic  boron  nitride 

pyrolytic  borcn  nitride 
phenolic  carbon 


phenolic  carbon 
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four  stations  are  presented  in  Figures  22  through  25.  It  is  interesting  to  ob¬ 
serve  that  the  region  of  greatest  chemical  interaction  appears  to  shrink  further 
and  further  away  from  the  outer  edge  of  the  boundary  layer  as  the  flow  progresses 
downstream.  This  is  consistent  with  the  knowledge  that  the  laminar  sublayer 
grows  at  a  slower  rate  than  the  turbulent  core  flow.  Another  interesting  feature 
of  these  species  profiles  is  in  Figure  25,  where  it  can  be  seen  that  a  signifi¬ 
cant  mole  fraction  of  boron  compounds  persists  in  the  boundary  layer  although 
the  flow  is  adjacent  to  a  phenolic  carbon  surface  at  that  station.  A  total  of 
60  species  were  considered  in  this  problem,  but  only  those  exhibiting  a  mole 
fraction  greater  than  0.001  somewhere  in  the  boundary  layer  are  shown  graphically. 
Finally,  Figure  26  presents  the  electron  concentration  profile  at  S  =  0.03231 
feet  over  the  BN  surface.  The  y-scale  is  left  unloqged  in  this  graph  in  order 
to  present  a  better  picture  of  the  actual  dimensions  of  the  regions  of  importance. 

3.  TRANSVERSE  CURVATURE  EFFECTS 

The  transverse  curvature  option  of  the  BLIMP  program  is  operational  for 
either  laminar  or  turbulent  reacting  flows.  The  simplest  of  these,  all  laminar 
flow,  was  chosen  for  a  sample  problem.  Air  was  assumed  to  be  flowing  over  a 
sharp  cone  at  the  following  conditions: 

11.27  atmospheres 
1650  Btu/lb 
6.38 

angle  =  9° 

Seven  nodes  were  chosen  through  the  boundary  layer,  and  the  problem  was  run  both 
with  and  without  consideration  of  transverse  curvature  (TVC) .  Some  of  the  re¬ 
sults  are  presented  in  Figures  27  through  29.  Figure  27  shows  the  TVC  effect  on 
drag  coefficient.  Inclusion  of  transverse  curvature  on  this  9°  ccne  increased 
Cf  by  as  much  as  25  percent.  A  comparison  of  Reynolds  number  on  momentum  thick¬ 
ness  is  presented  in  Figure  28.  Since  r/rQ  is  everywhere  greater  than  one,  one 
might  expect  from  the  definition  of  6  that  momentum  thickness  will  be  larger 
when  TVC  effects  are  considered.  The  steeper  slope  of  the  momentum  thickness 
curve  (de/ds)  for  this  zero  pressure  gradient  problem  is  also  apparent.  This 
increase  in  slope  is  expected  for  the  TVC  case  since  the  drag  coefficient  is 
larger.  Figure  29  shows  the  velocity  profile  comparison  at  S  =  0.025  feet. 
Transverse  curvature  effects  are  se  k  to  decrease  the  boundary  layer  thickness, 
thereby  increasing  wall  shear. 
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Figure  24.  Species  Concentration  Profiles  at  S  =  0.0321  Feet  for 
Sphere  Cone  Sample  Case 
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Figure  25.  Species  Concentration  Profiles  at  S  =  0.06273  Feet  for 
Sphere  Cone  Sample  Case 
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SECTION  V 

COUPLED  ABLATOR  BOUNDARY  LAYER  ENVIRONMENT  PROGRAM 

The  usual  ablating  heat  shield  or  nozzle  material  performance  analysis  re¬ 
quires  separate  determinations  of  the  inviscid  flow  field  over  the  body,  the 
boundary  layer  flow  near  the  surface,  and  the  surface  ind  in-depth  response  of 
the  ablating  material.  The  thought  has  no  doubt  crossed  the  mind  of  most  thermo¬ 
dynamics-fluid  dynamics  engineers  working  in  this  field  that  some  or  all  of  these 
analyses  should  be  combined  to  represent  the  coupled  physical  processes  more  cor¬ 
rectly  and  to  eliminate  some  of  the  wasted  effort  in  carrying  out  separate  analy¬ 
ses.  Indeed,  a  coupled  approach  is  desirable  since  the  material  response  affect; 
the  structure  of  the  boundary  layer,  and  the  boundary  layer  determines  the  <;ner< • 
and  mass  fluxes  at  the  material  surface  which  in  turn  control  the  heat  shield 
response.  Of  course,  any  change  in  body  shape  will  affect  the  inviscid  flow 
field.  The  totally  coupled  analysis  is,  unfortunately,  an  enormous  problem  if 
accurate  boundary  layer  and  material  response  calculations  aro  required,  hence 
the  present  separate  analysis  state  of  the  art.  Still,  it  is  possible  to  inves¬ 
tigate  coupled  charring  material,  multicomponent  boundary  layer  flow  problems 
within  che  context  of  certain  limiting  assumptions.  The  CABLE  program  incor¬ 
porates  subroutine  versions  of  the  BLIMP  program,  described  earlier  in  this  re¬ 
port,  and  the  Charring  Material  Ablation  (CMA)  program,  described  in  reference 
30,  to  accomplish  this  coupling.  An  earlier  version  of  the  CABLE  program  has 
been  described  elsewhere  (reference  31)  therefore  its  operation  will  only  be  sum¬ 
marized  in  Section  V.l.  Section  V.2  describes  recent  results  obtained  with  the 
program. 

1.  COUPLING  THE  CMA  AND  BLIMP  PROGRAMS 

The  CMA  program  is  an  implicit  finite  difference  one-dimensional  charring 
ablation  material  analysis  program  which  accounts  for  area  chan  ;e  due  to  mate¬ 
rial  curvature  in  a  general  fashion,  with  planar,  cylindrical,,  and  spherical 
geometries  as  special  cases.  Temperature  dependent  thermal  properties  aie  al¬ 
lowed,  and  the  user  may  specify  kinetically  controlled  pyrolysis  reactions.  Due 
to  its  one-dimensional  nature  however,  the  CMA  program  may  not  be  particularly 
accurate  in  regions  where  lateral  conduction  is  significant,  such  as  sharp  nose 
tips.  Al-,o,  since  local  static  pressure  is  input  into  the  BLIMP  program,  no 
account  is  taken  of  the  body  shape  change  on  the  pressure  distribution  i'u  tic 
coupled  analysis.  Within  these  limitations,  tne  CABLE  program  ■  ves  a  veiv  de¬ 
tailed,  accurate  picture  of  the  ablation-boundary  layer  interaction. 
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There  are  a  number  of  ways  in  which  the  boundary  layer  and  charring  abla¬ 
tion  programs  can  be  coupled.  In  the  procedure  which  has  been  adopted,  the  CMA 
program  is  effectively  the  controlling  program,  calling  BLIMP  whenever  necessary. 
Since  the  time  (0)  steps  in  a  charring  ablation  analysis  are  typically  quite 
small  (0.001-1.0  seconds),  yet  the  times  of  significant  change  in  a  boundary 
layer  flow  are  fairly  far  apart  (0.5-10.0  seconds),  it  would  be  uneconomical  to 
require  a  boundary  layer  solution  at  every  CMA  time  step.  Therefore  the  CABLE 
technique  is  to  input  mandatory  solution  times  at  which  complete  boundary  layer 
solutions  are  carried  out.  Between  these  mandatory  times  boundary  layer  param¬ 
eters  are  found  by  interpolation,  and  these  boundary  layer  parameters  are  pro¬ 
vided  to  the  CMA  program  as  boundary  conditions.  Another  computer  time  saving 
device  is  the  use  of  discrete  values  of  normalized  pyrolysis  gas  mass  flow  rates 

(m*)  and  normalized  char  mass  flow  rates  (m*) .  As  the  charring  ablation  program 

•  • 

proceeds  in  time,  a  new  combination  of  m*  and  m*  evolves  from  the  CMA  solution 
at  each  time  step.  It  would  not  be  possible  to  anticipate  and  calculate  the 
boundary  layer  flows  for  each  of  these  m*-m£  combinations  at  every  mandatory  so¬ 
lution  time,  therefore  a  grid  of  m*-m*  values  is  preselected  at  which  BLIMP  cal¬ 
culations  may  be  run.  The  BLIMP  program  is  automatically  run  at  several  of  the 
preselected  grid  points  in  the  assigned  injection,  surface  equilibrium  (or  kinet- 
ically  controlled  reaction)  mode,  and  the  required  boundary  layer  parameters  are 
found  from  a  three  way  interpolation  of  the  boundary  layer  solutions  in  the  m*, 

m*,  0  coordinate  array.  It  is  significant  that  not  all  BLIMP  solutions  in  the 

•  • 

m*,  m*,  0  array  need  be  run,  but  rather  only  those  which  are  needed  to  provide 
interpolation  values  for  the  current  CMA  solution.  Initially,  eight  BLIMP  solu¬ 
tions  are  required,  which  consist  of  two  values  of  m*  with  two  values  of  m*,  all 
at  the  first  two  mandatory  times.  The  CMA  solution  then  proceeds  until  the  solu¬ 
tion  "path"  exceeds  the  limits  of  the  "cube"  defined  by  the  first  eight  (m*,  m*, 

9  c 

0)  coordinate  points.  Additional  BLIMP  solutions  are  then  called  for  at  the  new 

required  grid  points  (m*,  m*,  0)  such  that  ordinary  interpolation  can  take  place, 

y  c 

whereupon  the  CMA  program  again  proceeds.  Thus,,  although  a  fairly  extensive 
•  • 

array  of  m*,  m*,  0  values  may  be  defined,  only  the  required  BLIMP  solutions  will 
9  c 

be  run.  This  procedure  will  be  explained  more  thoroughly  by  an  example  below. 

In  the  procedure  which  has  been  adopted,  the  transient  charring  ablation 
solution  is  effectively  the  controlling  program.  The  charring  ablation  solution 
at  a  given  station  proceeds  noniteratively ,  calling  the  boundary-layer  procedure 
as  needed  to  fill  in  the  surface  boundary  condition  matrix.  The  complete  time 
history  at  each  body  station  is  performed  prior  to  advancing  to  the  next  body 
station.  As  an  example  of  the  procedure,  consider  a  single  body  point  being  ana¬ 
lyzed  by  the  coupled  program  between  two  mandatory  solution  times  0^  and  G2>  The 
diagrams  below  indicate  the  projections  in  the  planes  0  =  0^  and  0  =  02  of  a 
hypothetical  history  of  m*  and  m*  as  generated  by  the  CMA  program  between  the  two 
mandatory  times. 
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0  1  2  3  4 

m* 


The  solution  at  time  0^  is  indicated  by  asterisks,  whereas  the  solution 
at  time  02  is  indicated  by  circles.  The  grid  values  for  m*  =  0,  1,  2,...  and 
m*  =  0,  1,  2,...  are  the  preselected  values  for  these  parameters  at  which  para¬ 
metric  boundary  layer  solutions  are  conducted  if  and  when  needed.  Based  on  the 
point  (*)  at  time  0^,  boundary- layer  solutions  are  generated  for  the  (m*,  m*) 
points  (1,1),  (1,2),  (2,1),  and  (2,2)  at  times  6-^  and  02-  Charring  ablation  so¬ 
lutions  can  be  obtained  for  times  01<0<02by  linear  interpolation  as  long  as 
m*  and  m*  stay  within  these  values.  Suppose  that  the  course  of  the  calculation 
between  times  0.^  and  02  xs  as  indicated  in  the  sketch.  Then,  additional  solu¬ 
tions  at  m*,  m*  of  (1,3)  and  (2,3),  then  (3,2)  and  (3,3)  and  finally  (2,4)  and 
y  c 

(3,4)  would  be  required,  each  at  both  times.  When  time  @2  (point  ©  )  is  ap¬ 
proached,  the  BLIMP  program  is  called  upon  for  a  solution  at  time  02  for  the 
exact  values  of  m*  and  m*  required  by  CMA.  This  boundary-layer  solution  is 
printed  out  and  that  information  needed  for  future  reference  (at  downstream  sta¬ 
tions)  is  saved  on  tape.  Solutions  are  then  performed  for  time  0^  for  the  cur¬ 
rent  bracketing  values  of  m*,  m*  (in  the  present  example,  values  of  (m*,  m*)  of 

9  c  g  c 

(2,3),  (3,3),  (2,4),  and  (3,4)).  These  boundary  layer  solutions  at  time  0^  are 
placed  over  those  for  0^  by  a  tape  flip-flop  since  the  latter  are  no  longer 
needed.  The  charring  ablation  solution  next  proceeds  from  time  ©2  to  time  03, 
calling  the  BLIMP  program  only  in  the  event  that  this  range  of  (m*,  m*)  is 
exceeded. 

The  above  described  procedure  is  over-simplified  to  some  extent  since  for 
many  materials  of  interest  in  certain  ablation  regions  the  use  of  m£  as  an  in¬ 
dependent  variable  is  a  poor  choice.  Consider  the  peculiarities  of  carbon  ablat¬ 
ing  with  air,  for  example,  as  illustrated  in  Figure  30,  which  shows  the  depen¬ 
dence  on  temperature  of  carbon  ablation  rates  in  air.  Ov»r  a  wide  range  of  tem¬ 
perature  m*  is  for  all  practical  purposes  independent  of  temperature;  in  this 
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range  the  only  useful  independent  variable  in  the  surface  thermochemical  solution 

process  is  temperature,  since  it  would  be  nearly  impossible  to  select  a  sequence 

of  m*  values  distributed  along  the  plateau.  At  high  temperatures,  however,  m* 
c  .  c 

becomes  strongly  dependent  on  T,  so  that  m*  is  the  preferred  independent  variable. 

c  , 

The  CABLE  program  logic  has  been  altered  recently  to  accept  either  m*  or  Tw  as 
an  independent  variable.  In  the  typical  problem  situation,  the  user  begins  his 
sequence  of  independent  variables  with  temperature  entries  starting  at  the  low¬ 
est  temperature  of  interest  and  extending  up  to  (as  a  minimum  objective)  the 
point  where  m*  begins  to  rise  above  the  plateau  value  for  the  lowest  pressure  to 
be  encountered.  It  would  be  conservative  and  preferable  to  add  temperatures  even 
beyond  the  "rise  point"  at  the  highest  pressure  in  case  these  points  should  be 
needed  during  the  solution.  After  these  temperature  values,  m*'s  are  added  to 


Figure  30.  Ablation  Rate  m*  Versus  Temperature  for  Carbon  in  Air 

span  the  expected  range  of  m*.  The  user  is  careful  to  pick  a  minimum  m*  slightly 
above  the  plateau  value.  An  m*  value  near  the  start  of  the  plateau  or  off  the 
plateau  in  the  low  temperature  region  is  extremely  undesirable  for  two  reasons: 

1.  The  lowest  assigned  m*  determines  the  break  between  assigned  temper- 

c  .  . 

ature  solutions  and  assigned  m*  solutions;  a  low  m*  will  in  effect 
"disqualify"  all  of  the  assigned  tempeiature  points  across  the  pla¬ 
teau  and  leave  a  large  "hole"  in  the  array  of  surface  thermochemical 
solutions . 

2.  For  equilibrium  calculations,  assigned  m*  solutions  at  m*  values  be¬ 
low  the  plateau  have  a  high  probability  of  nonconvergence. 

The  coupling  procedure  which  has  been  summarized  above  is  very  straight¬ 
forward  in  practice,  although  difficult  to  describe  in  words.  The  sample  prob¬ 
lem  presented  in  Section  v.2  should  clarify  some  of  the  more  difficult  points. 
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The  storage  requirements  for  this  coupling  procedure  are  surprisingly 
small.  In  the  first  place,  the  charring  ablation  solutions  are  noniterative  and 
the  complete  solution  for  all  times  at  any  one  station  is  accomplished  and  the 
results  printed  out  before  advancing  to  the  next  station.  Thus  no  historic  in¬ 
formation  relative  to  the  charring  ablation  solution  has  to  be  stored.  With  re¬ 
gard  to  the  boundary  layer,  only  two  mandatory  times  with  four  (m*,  m*)  combina¬ 
tions  at  each  of  these  times  need  to  be  considered  at  the  same  time.  The  only 
quantities  in  the  boundary  layer  which  need  to  be  dimensioned  for  the  full  time 
array  are  three  input  quantities  of  time,  total  pressure,  and  total  enthalpy. 

Edge  conditions  are  computed  around  the  body  at  the  time  of  the  stagnation-point 
calculation  since  the  necessary  integrations  are  performed  by  curve  fitting. 

This  necessitates  that  streamwise  dimension,  static  pressure,  edge  velocity,  edge 
density,  edge  viscosity,  edge  temperature,  body  curvature  (rQ) ,  transformed 
streamwise  dimension  (5),  pressure  gradient  parameter  ($) ,  and  the  flux  normaliz¬ 
ing  parameter  (a*)  be  dimensioned  for  the  number  of  streamwise  positions  (but 
not  for  time) .  About  300  numbers  must  be  stored  during  the  flip-flop  operation 
associated  with  the  two  times  which  are  being  considered  simultaneously,  whereas 
about  500  numbers  must  be  stored  on  tape  to  reenter  the  boundary  layer  at  the 
same  time  but  at  the  next  downstream  station  (used  for  first  guesses  and  for  cal¬ 
culation  of  nonsimilar  terms).  Thus,  both  permanent  machine  storage  requirements 
and  tape  storage  requirements  are  not  excessive  as  a  consequence  of  coupling. 


This  coupling  approach  has  the  important  feature  that  the  CMA  program  op¬ 
erates  very  nearly  as  it  does  when  used  in  conjunction  with  the  ACE  program  (see 
reference  32) .  In  the  CMA/ACE  approach,  complete  surface  tables  are  computed 
a  priori  and  independently  with  the  ACE  program  and  these  are  available  to  the 
CMA  solution.  In  the  coupled  approach,  these  surface  tables  are  initialized 
with  the  word  VOID.  When  the  CMA  program  encounters  this  word,  the  BLIMP  program 
is  called  to  supply  the  requisite  information  for  that  6,  m*  and  m*  (or  Ty) .  It 
is  thus  significant  that  the  CMA/ACE  approaches  have  been  used  extensively  and 
very  successfully  for  a  wide  variety  of  materials  and  environments.  Likewise, 
the  boundary- layer  calculations  are  performed  with  assigned  m*  and  m*  or  assigned 

m*  and  T  ,  together  with  the  requirement  of  surface  equilibrium  (with  possible 
g  w 

specified  rate-controlled  surface  reactions),  options  of  the  BLIMP  program  which 
also  have  been  exercised  extensively  with  success.  Furthermore,  this  replacement 
of  the  wall  mass  and  energy  balances  by  these  simple  assignment  statements  adds 
stability  to  the  boundary-layer  solution. 


2.  A  SAMPLE  PROBLEM  AND  RESULTS 

As  a  demonstration  of  the  CABLE  program,  a  coupled  transient  solution  was 
run  for  a  sphere-cone  reentry  body  with  a  nose  radius  of  0.5  inches  and  a  cone 
half-angle  of  eight  degrees.  A  single  body  point  at  a  surface  running  length  of 
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10  inches  (based  on  the  original  body  shape)  was  analyzed  although  there  is  no 
reason  other  than  computer  time  limitations  why  more  body  stations  could  not 
have  been  specified.  Laminar  flow  was  assumed  for  the  boundary  layer  to  hold 
the  required  number  of  nodes  to  a  minimum.  The  heat  shield  material  was  one- 
inch  thick  phenolic  carbon  with  properties  as  described  in  reference  28  and  was 
insulated  at  the  back  face.  The  vehicle  trajectory  is  shown  in  Figure  31.  Total 
flight  time  from  300,000  feet  is  27  seconds  with  an  impact  velocity  of  17,000  ft/ 
sec.  Pressure  ratio  P/PQ  at  the  body  station  of  interest  was  assumed  to  be 
0.0262  for  the  entire  27  seconds.  Mandatory  BLIMP  solution  times  were  0.0,  7.0, 
14.0,  18.0,  21.0,  23.0,  25.0,  and  27.0  seconds.  A  total  of  144  BLIMP  solutions 
in  the  assigned  injection,  surface  equilibrium  mode  exclusive  of  these  mandatory 
solutions  were  run  to  provide  the  CMA  program  with  boundary  conditions. 

Typical  CABLE  results  are  shown  in  Figures  32-40.  Surface  temperature 
history  is  shown  in  Figure  32,  where  it  is  seen  that  a  maximum  temperature  of 
4480°R  was  reached  at  the  end  of  the  flight.  Figure  33  illustrates  the  pyroly¬ 
sis  gas  flow  rate  history  at  the  surface.  Little  or  no  pyrolysis  occurs  during 
the  first  ten  seconds  of  flight,  with  a  peak  outgassing  rate  of  approximately 
0.028  lbs/ft2sec  reached  near  the  end  of  the  flight.  Figure  34  describes  the 
progression  of  the  pyrolysis  and  char  fronts  into  the  heat  shield  material.  The 
pyrolysis  front  is  arbitrarily  defined  as  the  point  where  the  local  material  den¬ 
sity  is  equal  to  the  char  density  plus  98  percent  of  the  difference  between  the 
virgin  and  char  densities.  The  char  front  is  similarly  defined  as  that  point 
where  the  local  material  density  is  equal  to  the  char  density  plus  2  percent  of 
the  difference  between  the  virgin  and  char  densities.  The  pyrolysis  front  pene¬ 
trates  0.208  inches  of  the  heat  shield  during  this  trajectory,  whereas  the  fully 
charred  material  reaches  a  depth  of  only  0.067  inches  by  the  end  of  the  flight. 

No  surface  recession  occurs  at  this  body  station. 

Figures  35-37  show  profiles  for  some  of  the  thermodynamic  variables  of 
interest  for  this  particular  body  station  at  t  =  21  seconds.  The  subsurface 
density  profile  of.  Figure  35  indicates  the  extent  of  the  pyrolysis  region.  Fig¬ 
ure  36  contains  che  velocity  profile  for  the  boundary  layer  flow,  while  Figure  37 
presents  temperature  profile  information.  The  coupling  of  the  boundary  layer  to 
the  subsurface  material  is  particularly  apparent  in  the  surface  temperature  con¬ 
tinuity  illustrated  in  Figure  37.  Figures  38-40  contain  the  species  mole  fraction 
profiles  through  the  boundary  layer  at  three  times  of  interest.  At  t  =  0.0 
seconds,  the  cold  wall  forces  a  recombination  of  N  atoms  to  0  atoms  to  0 2, 
etc.,  giving  the  unusual  profiles  of  Figure  38.  At  t  =  21  seconds  (Figure  39) 
significant  pyrolysis  has  occurred  with  the  resulting  gases  being  injected  into 
the  boundary  layer  flow.  Hydrogen,  carbon,  and  species  containing  these  atoms 
are  evident  in  the  boundary  layer  gas.  At  t  =  27  seconds  (Figure  40)  the  boundary 
layer  species  are  much  the  same,  however  the  lower  edge  temperature  caused  by  the 
slowing  down  of  the  vehicle  has  allowed  many  more  trace  species  to  be  formed. 
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Figure  31.  Trajectory  for  CABLE  Sample  Problem 
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Figure  32.  Surfa< 
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Figure  33.  Pyrolysis  Gas  Flow  Rate  History  at 
S  =  10.0  Inches 
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Figure  39.  Mole  Fraction  Profiles  at  S  =  10.0  Inches, 
t  =  21.0  Seconds 
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The  above-described  CABLE  run,  comprised  of  a  total  of  152  BLIMP  solutions 
plus  the  in-depth  charring  material  solution,  required  approximately  30  minutes 
on  a  UNIVAC  1108  computer.  Solutions  at  succeeding  body  stations  should  proceed 
somewhat  faster  due  to  the  better  initial  guesses  provided  by  the  upstream  solu¬ 
tion,  however  if  turbulent  flow  situations  are  expected,  about  twice  as  many 
boundary  layer  nodes  would  be  needed.  Thus,  it  is  estimated  that  a  complete  re¬ 
entry  vehicle  body  solution  with  the  CABLE  program,  including  29  body  stations  in 
a  27  second  trajectory,  would  require  from  10  to  20  hours  of  UNIVAC  1108  computer 
time. 
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SECTION  VI 

CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER  WORK 

The  BLIMP  boundary  layer  code  has  been  successfully  extended  to  include 
turbulent  flows  and  flows  with  transverse  curvature.  The  basic  analytical  tech¬ 
nique  offers  the  advantages  of  completely  general  multicomponent  chemistry,  pre¬ 
servation  of  nonsimilar  terms  in  the  equations  of  motion,  and  intimate  coupling 
to  the  surface  material  behavior.  The  coupling  occurs  through  requirements  of 
surface  equilibrium  and/or  a  surface  energy  balance,  or  a  complete  coupling  with 
an  in-depth  conduction  and  charring  ablation  code. 

There  are  several  areas  where  the  code  and  the  analytical  model  which  it 
includes  can  be  improved  or  extended.  One  of  these  is  in  the  general  area  of 
modeling  turbulence  in  the  equations  of  motion.  The  purpose  of  the  work  pre¬ 
sented  in  this  report  was  to  incorporate  a  turbulent  model  into  the  BLIMP  pro¬ 
gram.  The  model  that  was  used  is  perhaps  the  best  available  for  general  hyper¬ 
sonic  flows  with  arbitrary  species  injection,  however  an  extensive  investigation 
into  the  generality  of  the  model  in  all  flow  situations  is  called  for.  Experi¬ 
mental  data  comparisons  at  high  Mach  numbers  and  possible  model  changes  to  match 
these  data  would  be  most  valuable. 

Two  smaller  changes  in  the  program  would  make  it  more  accurate  and  conven¬ 
ient  to  use  on  some  problems.  The  first  of  these  involves  the  aH  or  coordinate 
stretching  parameter.  On  very  long  bodies,  the  boundary  layer  thickness  varies 
over  two  or  three  orders  of  magnitude.  The  aH  parameter  forces  the  boundary 
layer  thickness  to  remain  constant  in  the  solution  plane,  thereby  eliminating 
the  need  for  a  large  number  of  nodes  which  are  unused  near  the  stagnation  point 
or  leading  edge.  For  turbulent  flows,  the  laminar  sublayer  grows  at  a  much 
slower  rate  than  the  total  boundary  layer.  As  the  normal  coordinate  is  stretched 
in  the  solution  plane,  the  nodes  nearest  the  wall  eventually  can  be  pulled  out¬ 
side  the  sublayer  altogether,  resulting  in  poor  convergence  of  the  numerical 
technique  and  inaccurate  solutions.  A  different  stretching  technique  scaled  to 
the  sublayer  as  well  as  the  boundary  layer  edge  must  be  invented  to  avoid  this 
problem.  The  second  change  involves  nonisentropic  expansions.  It  is  currently 
possible  to  specify  entropy  and  pressure  in  order  to  fix  the  edge  state  of  the 
boundary  layer  gas  for  "entropy  layer"  flows.  This  technique  introduces  some 
inaccuracy  however  since  an  entropy  gradient  in  the  streamwise  direction  also  re¬ 
quires  a  velocity  gradient  in  the  normal  direction  (nonzero  vorticity) .  The 
analysis  should  be  modified  to  include  this  edge  velocity  gradient  for  nonisen¬ 
tropic  expansions. 
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The  final  recommendation  is  in  the  general  category  of  electron  collision 
models.  At  present,  a  very  simple  clean  air  collision  frequency  model  is  in¬ 
cluded  in  the  program,  and  electron  distributions  are  calculated  assuming  equi¬ 
librium  chemistry.  More  accurate  clean  air  collision  frequency  models  are  avail¬ 
able  and  could  be  incorporated  into  the  program.  However,  based  on  the  level  of 
sophistication  of  the  rest  of  the  boundary  layer  analysis,  it  seems  appropriate 
to  attempt  to  take  into  account  the  effect  of  the  ablation  products  on  the  elec¬ 
tron  collision  frequencies.  Also,  charge  separation  effects  may  alter  the  equi¬ 
librium  electron  distribution  significantly.  It  is  recommended  that  a  study  be 
undertaken  to  establish  the  importance  of  charge  separation  and  ablation  products 
on  the  basic  parameters  of  interest  in  communications  and  that  these  effects  be 
included  in  the  BLIMP  code  if  necessary. 
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APPENDIX  I 

COORDINATE  TRANSFORMATIONS 

A  convenient  coordinate  transformation  for  use  in  boundary  layer  problems 
with  transverse  curvature,  and  the  one  used  in  the  present  analysis,  can  be 
arrived  at  by  nondimensionalizing  and  simplifying  the  continuity  and  momentum 
equations  in  a  relatively  straightforward  manner.  The  momentum  equation  for  a 
turbulent  flow  with  transverse  curvature  is 


K  3u  ,  K  3ll  3  r  K  , 
pur  +  pvr  gy  =  gy  [pr  (v  +  c 


M;  3yJ 


(146) 


while  the  requirement  of  continuity  yields 


3purK  ,  3pvrK  « 
3s +  3y  0 


(147) 


The  two  equations  can  be  conveniently  combined  by  the  definition  of  a  stream 
function  f,  with  an  arbitrary  wall  value,  f  : 


ry 

_  /  pur 
L  prurro 


(148) 


The  stream  function  is  made  dimensionless  by  the  introduction  of  reference  condi¬ 
tions  p  and  u  and  by  the  v-dimension  scaling  parameter  6  (x) .  The  reference 
condition  (  )  will  be  taken  as  the  isentropic  boundary  layer  edge  condition, 
while  the  scaling  length  <5  remains  to  be  defined.  Solving  the  continuity  equa¬ 
tion  for  pvrK,  it  is  easily  shown  that 


pvr<  =  pwvwro  -  h 


rjprur6(f  -  fw)_ 


(349) 


Using  this  result  in  the  momentum  equation  yields 


5?  *  !'.V!  -  h  co-Iur5 


<*  -  v  S? 


3  r  < 

-  n  r 
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.  ,  „  >  3u  k  3P 
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Introducing 


The  momentum  equation  becomes 


pGrKur  [a  oT  +  Ur  If  ‘  Ur  6  H  “  j  +  jpwVwro 

-  Is  [roprur6  (f  *  V]  +  6  H  ~  [roprur6  (f 
1  3  f"  k  Ur  .  ,  .  r  v  3ul  K  dP 

■  * «  L  T  M  »J  E 


,11  fUr  3u"! 


(152) 


A  new  normal  coordinate  is  now  defined  which  allows  a  simple  relationship  be¬ 
tween  the  dimensionless  stream  function  f  and  dimensionless  velocity  u.  Noting 


3f  pur 


3y  Prurrc 


then  if  we  define 


dn 

H  prro 


the  relation  between  f  and  u  becomes 


3f  ,, 

37)  =  f  °KU 


The  au  is  a  stretching  parameter  on  the  normal  coordinate  ri  such  that  the  bound- 
ary  layer  is  contained  within  a  constant  n  range,  0  to  ne<jge  (see  Section  II.  3). 
Using  the  new  n  coordinate  and  expressing  u  in  terms  of  f',  the  momentum  equa¬ 
tion  takes  the  form 
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An  additional  coordinate  transformation  will  bring  the  equation  to  a  more 
pact  form.  Defining 

d£  =  Qds 

where  Q(s)  is  as  yet  undefined,  the  streamwise  derivatives  become 


d  (  )  Q  d  (  ) 
ds  E,  d  In  £ 


d  ln  ur  q  d  ln  ur 
ds  £  d  in  £ 


QL 


which  leads  to 


/  P,.\ 

!(f'2  “  aH  T) 


+  2f 


,  3f ' 


3  ln  £ 


2f 


d  ln  a„ 

i  2  _ H 

d  ln  X 


-  f" 


3f 


3  ln  X 


+  2f 


d  ln(rj  u  i) 
o  r  r 

2prS 

t(C  +  eM)f" 

d  ln  C 

Q52Prur 

aH 

Defining  Q  and  6  properly  will  yield  unity  for  two  groups  of  terms  in  the 
version  of  the  equation.  Choosing  the  group  in  the  left  hand  side  first 

d  ln(/p  uj) 


which  requires 


The  second  group  of  terms  is 


6  = 


roP  U 

o'r  r 


2Mr" 

Q6 1  p  a 
r  r 


=  1 


which  requires 


O  =  o  u  u  r2  ** 
v  -  ur  rMr  0 


With  tnese  changes,  the  erne '  turn  equation  becomes 


com- 

(163) 

(164) 

(165) 

(166) 
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(167) 

(.16  8) 
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rl 


(171) 


The  final  coordinate  transformation  is 


£  =  /  Pup  r2K 

/  ^r  ryr  o 

Jo 


ds 


(172) 


n .  fy 

Jo  prro 


dy 


(173) 


while  the  other  definitions  take  the  form 


purKdy 


(174) 


s 

f  =  -  —  /  p  v  rKds 

w  m  J  w  w  ° 


(175) 


-103- 


AFWL-TR-69-106 


APPENDIX  II 

DEVELOPMENT  OF  A  RECURSION  FORMULA 
FOR  MIXING  LENGTH 


The  mixing  length  differential  equation  has  been  shown  to  be  of  the  form 


where 


The  general  solution  is 


=  (KV  “  l)p 


P(n)  =  -^4 - 


l  =  KaH(r)  -  L) 


(178) 


where 


•n 


L(n)  =  -- 


f  *  ° 


/  '  Pdn1 

*b 

e 


In  the  definition  of  L  at  the  i'th  rode,  the  integrals  from  0  to  can  be 
broken  into  two  parts,  0  to  and  to  to  give 


(\  ^ 

Li-1  A  7n i-1 


Pdn' 

i  ■?  —  1 


rni 

/  Pdn 


(180) 


X  1  Pdn 
*Vl 


Assuming  a  linear  variation  of  P(n)  over  the  interval  ni_1  to  r>i,  the  first 


term  of  the  Lj,  expression  becomes 
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This  last  integral  term  can  be  broken  into  two  pieces  to  give 


n  fni  Pdn ' 

rn-  J 

J  e  Vi  dn 
ni-l 


-d2 
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hi 


1  Pdn 
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=  A[Dw(d)  -  BDw(c)] 


(188) 


The  Dawson  Integral,  Dw(  )*,  can  be  evaluated  from  tables  (Ref.  23)  or,  in  the 
case  of  the  present  analysis,  by  a  series  method.  The  quantity  is  given  by 


L. 

l 


(189) 


which  is  a  recursion  formula  for  in  terms  of  the  value  for  L  at  the  previous 
node  and  the  local  values  of  P(n). 
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